//
//-----------------------------------------------------------------------
/// @copyright
///            (c) Copyright 2008 by GATS, Inc.,
///     11864 Canon Blvd, Suite 101, Newport News VA 23606
///
///  All Rights Reserved. No part of this software or publication may be
///  reproduced, stored in a retrieval system, or transmitted, in any form
///  or by any means, electronic, mechanical, photocopying, recording, or
///  otherwise without the prior written permission of GATS, Inc.
///
//-----------------------------------------------------------------------
///
/// @file	SunSensorCalc.cpp
///
/// @author	John Burton	<j.c.burton@gats-inc.com>
///
/// @date  	Thu Jan 17 14:08:25 2008
///
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//
#include <map>
#include <string>
#include "SunSensorCalc.h"
#include "EdgeModel.h"
#include "EdgeLocation.h"
#include "DriftCorrection.h"
#include "GATS_Utilities.hpp"
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//

static const double CSUM_INT    = 347.00; // Center Sum integration time (us)
static const double XROI_INT    = 467.67; // X Edge ROI integration time (us)
static const double YROI_INT    = 456.83; // Y Edge ROI integration time (us)
static const double MISSING     = -1e24;
static const double EDGE_FACTOR = 0.198847;
static const double FIXED2FLOAT = 32.0;
static const double ARCMIN_TO_RADIANS = M_PI / (double)10800.0;
//static double det_el = 636.15625;  // now in units of pixels
static double det_el = 639.039083;  // now in units of pixels
static double det_az = 518.46875;   // now in units of pixels

//
//-----------------------------------------------------------------------
// Class Methods:
//-----------------------------------------------------------------------
//


void SunSensorCalc::initNames(void)
{
  SunSensorTimeGrid_.setName("SolarImageTimes");
  ScienceTimeGrid_.setName("SS_ScienceTimeGrid");
  Track_Low_X_.setName("SolarTrackXLow");
  Track_High_X_.setName("SolarTrackXHigh");
  Track_Low_Y_.setName("SolarTrackYLow");
  Track_High_Y_.setName("SolarTrackYHigh");
  
  SolarHighEl_.setName("SolarHighElevationEdge");
  SolarLowEl_.setName("SolarLowElevationEdge");
  SolarHighAz_.setName("SolarHighAzimuthEdge");
  SolarLowAz_.setName("SolarLowAzimuthEdge");
  SolarExtent_.setName("SolarExtent");
  DetectorAz_.setName("DetectorLockdown_AZ");
  DetectorEl_.setName("DetectorLockdown_EL");
  RecordType_.setName("SolarTrackRecordType");
  CenterOffsets_.setName("SolarCenterSumsOffsets");

  NormalizedSolarExtent_.setName("NormalizedSolarExtent");
  
  SumsData_.setName("SolarSums");
  LocsVect_.setName("SolarSumsLocations");
  CenterSumsVect_.setName("SolarCenterSums");
  CenterRowsVect_.setName("SolarCenterRows");
  CenterAvesVect_.setName("CenterAveragePixelValues");
  AvePixVect_.setName("AveragePixelValues");
  CenterOffsetsVect_.setName("SolarCenterSumsOffsetsVect");
  
  QualityOfFit_.setName("SunSensorQualityOfFit");

  zImpact_.setName("ImpactAltitude");
}



EventVar<double> SunSensorCalc::cleanupTransmission(EventVar<double> &trans, EventVar<double> &ia)
{
  EventVar<double> tmptrans(trans);
  for(int i=0; i<tmptrans.size(); i++)
  {
    if(ia[i] < 0.35)
      tmptrans[i] = 1.0;
    else
      if((ia[i] >= 0.35)&&(ia[i] < 0.37))
      {
	double y = 18.5 - 50.0 * ia[i];
	tmptrans[i] = y + (1-y)*tmptrans[i];
      }
  }

  tmptrans = tmptrans.Smooth(21);
  tmptrans = tmptrans.Smooth(21);
  return tmptrans;
}


double SunSensorCalc::calcMaxIntensity(void)
{
  double tmp = 0.0;
  double maxint = 0.0;
  EventVar<double> avepix;
  valarray<bool> good  = ((SunSensorTimeGrid_ >= 60.0) &&
			  (SunSensorTimeGrid_ <= 80.0));

  for(int i=0; i<AvePixVect_.size(); i++)
  {
    avepix = AvePixVect_[i][good];
    tmp    = avepix.max();
    maxint = (maxint > tmp) ? maxint : tmp;
  }

  return maxint;
}


SunSensorCalc::SunSensorCalc(Event& L0, Event& L1, Event& Tmp, ConfigFile& cf, std::string Section)
{

  cf_      = cf;
  section_ = Section;
  srssflag_ = L0.getSRSSFlag();
  EventVarVect<int> intSums;
  EventVar<double>tempExtent;
  EventVarVect<double> tmpTrans;
  EventVarVect<double> tmpImpact;

  Tmp.getEventVar("SolarImageTimes",          SunSensorTimeGrid_);
  Tmp.getEventVar("SS_ScienceTimeGrid",       ScienceTimeGrid_);
  Tmp.getEventVar("SS_CompositeTransmission", tmpTrans);
  Tmp.getEventVar("SS_CompositeImpactAngles", tmpImpact);
  Tmp.getEventVar("SS_TopEdgeImpactAngle",    TopEdgeImpactAngle_);
  Tmp.getEventVar("SS_BottomEdgeImpactAngle", BottomEdgeImpactAngle_);
  Tmp.getEventVar("AveragePixelValues",       AvePixVect_);
  Tmp.getEventVar("SolarSumsLocations",       LocsVect_);
  Tmp.getEventVar("SolarTrackXLow",           Track_Low_X_);
  Tmp.getEventVar("SolarTrackXHigh",          Track_High_X_);
  Tmp.getEventVar("SolarTrackYLow",           Track_Low_Y_);
  Tmp.getEventVar("SolarTrackYHigh",          Track_High_Y_);
  Tmp.getEventVar("SolarCenterSums",          CenterSumsVect_);
  Tmp.getEventVar("CenterAveragePixelValues", CenterAvesVect_);

  L1.getEventVar("zImpact",        zImpact_);
  L1.getEventVar("L1_SolarExtent", tempExtent);

  TrueSolarExtent_ = tempExtent[0];

  Tmp.getEventVar("SolarExtent", tempExtent);

  tempExtent.setName("FromFPACorrection");

  SolarExtentCompare_.addEventVar(tempExtent);

  valarray<bool> valid  = (tmpTrans[0] > MISSING);
  
  Transmission_ = tmpTrans[0][valid];
  ImpactAngle_  = tmpImpact[0][valid];
  Transmission_ = cleanupTransmission(Transmission_,ImpactAngle_);
//  Transmission_.Plot("Transmission Vs. Impact Angle","",ImpactAngle_);

  eventNumber_    = L0.getEventNumber();
  MaxAveIntensity = calcMaxIntensity();

  if(cf_.GetFlag(section_,"DebugPrint"))
  {
    std::cerr << "SunSensorTimeGrid_.size()     = " << SunSensorTimeGrid_.size() << "\n";
    std::cerr << "ScienceTimeGrid_.size()       = " << ScienceTimeGrid_.size() << "\n";
    std::cerr << "Transmission_.size()          = " << Transmission_.size() << "\n";
    std::cerr << "ImpactAngle_.size()           = " << ImpactAngle_.size() << "\n";
    std::cerr << "zImpact_.size()               = " << zImpact_.size() << "\n";
    std::cerr << "TopEdgeImpactAngle_.size()    = " << TopEdgeImpactAngle_.size() << "\n";
    std::cerr << "BottomEdgeImpactAngle_.size() = " << BottomEdgeImpactAngle_.size() << "\n";
    std::cerr << "AvePixVect_.size()            = " << AvePixVect_.size() << "\n";
    std::cerr << "LocsVect_.size()              = " << LocsVect_.size() << "\n";
    std::cerr << "AvePixVect_[0].size()            = " << AvePixVect_[0].size() << "\n";
    std::cerr << "LocsVect_[0].size()              = " << LocsVect_[0].size() << "\n";
    std::cerr << "Track_Low_Y_.size()           = " << Track_Low_Y_.size() << "\n";
    std::cerr << "Track_High_Y_.size()          = " << Track_High_Y_.size() << "\n";
    std::cerr << "CenterSumsVect_.size()        = " << CenterSumsVect_.size() << "\n";
    std::cerr << "CenterAvesVect_.size()        = " << CenterAvesVect_.size() << "\n";
  }
  
  if(cf_.GetFlag(section_,"WriteScienceEdgeImpactAngle"))
  {
    std::ofstream outfile;
    std::string fname = "ScienceEdgeImpactAngle" + GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
    outfile.open(fname.c_str(),std::ofstream::out);

    int len = ScienceTimeGrid_.size();
    for(int i=0;i<len; i++)
    {
      outfile << i << ", " << ScienceTimeGrid_[i] << ", ";
      outfile << TopEdgeImpactAngle_[i] << ", " << BottomEdgeImpactAngle_[i];
      outfile << "\n";
    }
    outfile.close();
  }

  TopEdgeIntensity_      = calcEdgeIntensity(TopEdgeImpactAngle_);
  BottomEdgeIntensity_   = calcEdgeIntensity(BottomEdgeImpactAngle_);
  TopEdgeImpactAngle_    = TopEdgeImpactAngle_.VInterpol(ScienceTimeGrid_,SunSensorTimeGrid_);
  BottomEdgeImpactAngle_ = BottomEdgeImpactAngle_.VInterpol(ScienceTimeGrid_,SunSensorTimeGrid_);
  zImpact_               = zImpact_.VInterpol(ScienceTimeGrid_,SunSensorTimeGrid_);
  
  EventVar<double> diff;
  diff = BottomEdgeImpactAngle_ - TopEdgeImpactAngle_;
//  diff.Plot("Difference in Impact Angle","Bottom Edge - Top Edge");
  
  if(cf_.GetFlag(section_,"WriteSunSensorEdgeImpactAngle"))
  {
    std::ofstream outfile;
    std::string fname = "SunSensorEdgeImpactAngle" + GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
    outfile.open(fname.c_str(),std::ofstream::out);
    EventVar<double> timegrid;
    
    int len = SunSensorTimeGrid_.size();
    for(int i=0;i<len; i++)
    {
      outfile << i << ", " << SunSensorTimeGrid_[i] << ", ";
      outfile << TopEdgeImpactAngle_[i] << ", " << BottomEdgeImpactAngle_[i];
      outfile << "\n";
    }
    outfile.close();
  }

  if(cf_.GetFlag(section_,"WriteTransmission"))
  {
    std::ofstream outfile;
    std::string fname = "Transmission" + GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
    outfile.open(fname.c_str(),std::ofstream::out);
    int len = Transmission_.size();
    for(int i=0;i<len; i++)
    {
      outfile << i << ", ";
      outfile << Transmission_[i] << ", " << ImpactAngle_[i] << ", ";
      outfile << "\n";
    }
    outfile.close();
  }
}


void SunSensorCalc::writeExoAtmosphericDataXY(void)
{
  const int nelem = 7;
  int top[nelem]    = {0,1,2,3,4,5,6};
  int bottom[nelem] = {7,8,9,10,11,12,13};
  int center[nelem] = {24,25,26,27,28,29,30};
  int len = SunSensorTimeGrid_.size();
  int mid = len / 2;
  int first = mid - 200;
  int last  = mid + 200;
  int l;
  double edge_max_int =  MaxAveIntensity * XROI_INT / CSUM_INT; 

  std::string fname;
  std::string base = GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
  std::ofstream outfile;

  fname = "TopEdge_" + base;
  outfile.open(fname.c_str(),std::ofstream::out);
  for(int i=first; i<=last; i++)
  {
    for(int j=0; j<nelem; j++)
    {
      l = top[j];
      outfile <<  LocsVect_[l][i] << ", " << (AvePixVect_[l][i] / edge_max_int) << "\n";
    }
  }
  outfile.close();

  fname = "BottomEdge_" + base;
  outfile.open(fname.c_str(),std::ofstream::out);
  for(int i=first; i<=last; i++)
  {
    for(int j=0; j<nelem; j++)
    {
      l = bottom[j];
      outfile <<  LocsVect_[l][i] << ", " << (AvePixVect_[l][i] / edge_max_int) << "\n";
    }
  }
  outfile.close();

  fname = "Center_" + base;
  outfile.open(fname.c_str(),std::ofstream::out);
  for(int i=first; i<=last; i++)
  {
    for(int j=0; j<nelem; j++)
    {
      l = center[j];
      outfile <<  LocsVect_[l][i] << ", " << (AvePixVect_[l][i] / MaxAveIntensity) << "\n";
    }
  }
  outfile.close();
}



void SunSensorCalc::fillTrackingValues(void)
{
  const int nelem  = 14;
  int nsamp = AvePixVect_[0].size();
  int nrows = 1024;
  int l;
//  double edge_max_int =  MaxAveIntensity * XROI_INT / CSUM_INT; 
  double edge_max_int =  MaxAveIntensity; 
  EventVar<double> sample(nrows);
  EventVar<double> XAxis(0.0,(double)(nrows-1),1.0);

if(!cf_.GetFlag(section_,"FullElevationScanCalibrationEvent"))
{
  for(int i=0; i<nsamp; i++)
  {
    sample.setValueToMissing();
    for(int j=0; j<nelem; j++)
    {
      l = LocsVect_[j][i];
      sample[l] = AvePixVect_[j][i] / edge_max_int;
    }
    Track_Low_X_[i]  = findLowXEdge(TopEdgeIntensity_[i],sample,XAxis);
    Track_High_X_[i] = findHighXEdge(BottomEdgeIntensity_[i],sample,XAxis);
  }
}
}



void SunSensorCalc::fillFPATimeSeries()
{
  std::ofstream outfile;
  std::string fname;
//  const int nelem = 21;
  const int nelem = 14;
  int lt, lb, l;
//  int elems[nelem] = {0,1,2,3,4,5,6,24,25,26,27,28,29,30,7,8,9,10,11,12,13};
  int elems[nelem] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13};
  int nsamp = AvePixVect_[0].size();
  EventVar<double> Vals(nelem);
  EventVar<double> Locs(nelem);
  EventVar<double> EdgeData;
  QualityOfFit_.resize(nsamp);
//  EdgeModel SolarEdge(cf_,section_);
  double maxia        = ImpactAngle_.max();

  if(cf_.GetFlag(section_,"WriteEdgeRegressionParams"))
  {
    fname = "./EdgeRegressionParams_" + GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
    outfile.open(fname.c_str(),std::ofstream::out);
  }
//  SolarEdge.setTransmissionModel(ImpactAngle_,Transmission_);
  for(int i=0; i<nsamp; i++)
  {
    for(int j=0; j<nelem; j++)
    {
      l = elems[j];
      Locs[j]    = LocsVect_[l][i];
      Vals[j]    = AvePixVect_[l][i] / MaxAveIntensity;
    }
    EdgeModel SolarEdge(Vals,Locs,Track_Low_X_[i],Track_High_X_[i]);
//    EdgeData    = SolarEdge.findEdge();
//    if(BottomEdgeImpactAngle_[i] < maxia)
    SolarEdge.setTransmissionModel(ImpactAngle_,Transmission_);
    SolarEdge.attenuateEdgeModel(TopEdgeImpactAngle_[i],BottomEdgeImpactAngle_[i]);
    EdgeData = SolarEdge.findEdge(Vals,Locs,Track_Low_X_[i],Track_High_X_[i]);
    if(EdgeData[2] < 0.95)
    {
//      std::cerr << i << " - repeating SolarEdge.findEdge - fit quality = " << EdgeData[2] << "\n";
      SolarEdge.attenuateEdgeModel();
      EdgeData = SolarEdge.findEdge(Vals,Locs,Track_Low_X_[i],Track_High_X_[i]);
    }

//    if(eventNumber_ > 2564)
    if(!cf_.GetFlag(section_,"FullElevationScanCalibrationEvent") 
        && eventNumber_ > 2564)

    {
      Track_Low_X_[i]  = EdgeData[0];
      Track_High_X_[i] = EdgeData[1];
      QualityOfFit_[i] = EdgeData[2];
    }

    if(cf_.GetFlag(section_,"WriteEdgeRegressionParams"))
    {
      outfile << i;
      for(int j=3; j<10; j++)
	outfile << ", " << EdgeData[j];
      outfile << "\n";
    }
    if(cf_.GetFlag(section_,"WriteSampleModelData"))
    {
      std::ofstream fitfile;
      std::string name = "./SampleData/Model_" + GATS_Utilities::ConvertToString(i) + ".dat";
      fitfile.open(name.c_str(),std::ofstream::out);
      for(int j=0; j<nelem; j++)
	fitfile << Locs[j] << ", " << Vals[j] << ", " << EdgeData[j+10] << "\n";
      fitfile.close();
    }
  }
  if(cf_.GetFlag(section_,"WriteEdgeRegressionParams"))
    outfile.close();
  if(cf_.GetFlag(section_,"PlotQualityOfFit"))
    QualityOfFit_.Plot("Quality of Fit","");
  if(cf_.GetFlag(section_,"PlotEdgeLocationFPA"))
  {
    Track_Low_X_.Plot("Low X Edge","");
    Track_High_X_.Plot("High X Edge","");
    EventVar<double>temp;
    temp = Track_High_X_ - Track_Low_X_;
    temp.Plot("Diff","");
  }
}


  
EventVar<double> SunSensorCalc::calcEdgeIntensity(EventVar<double>& edge)
{
  EventVar<double> tmpIA;
  EventVar<double> tmpTrans;

  tmpTrans = Transmission_.VInterpol(ImpactAngle_,edge,1);
  tmpTrans = tmpTrans.VInterpol(ScienceTimeGrid_,SunSensorTimeGrid_);
  tmpTrans = tmpTrans * EDGE_FACTOR;

  return tmpTrans;
}




void SunSensorCalc::transformCoordinates(Event &SD)
{

  std::string fname;

  std::string Xname  = "X_Pixel_Spacing";
  std::string Xfile  = "X_MappingFile";

  std::string Yname  = "Y_Pixel_Spacing";
  std::string Yfile  = "Y_MappingFile";

  std::string det_az_name;
  std::string det_el_name;

  if(srssflag_ == 'r')
  {
    det_az_name = "Detector_lockdown_az_SR";
    det_el_name = "Detector_lockdown_el_SR";
  }
  else if(srssflag_ == 's')
  {
    det_az_name = "Detector_lockdown_az_SS";
    det_el_name = "Detector_lockdown_el_SS";
  }
  else
  {
    det_az_name = "Detector_lockdown_az";
    det_el_name = "Detector_lockdown_el";
  }
  EventVar<double> X_Mapping(Xname,511.0,-512.0,-1.0);
  EventVar<double> Y_Mapping(Yname,-511.0,512.0,1.0);

  if(SD.EventVarExists(Xname))
    SD.getEventVar(Xname,X_Mapping);
  else
  {
    if(cf_.ValidEntry(section_,Xfile))
    {
      fname  = cf_.GetStr(section_,Xfile);
      std::cerr << "Reading " << fname << std::endl;
      X_Mapping.readVar(fname);
    }
    X_Mapping = X_Mapping * ARCMIN_TO_RADIANS;
    SD.addEventVar(X_Mapping);
  }

  if(SD.EventVarExists(Yname))
    SD.getEventVar(Yname,Y_Mapping);
  else
  {
    if(cf_.ValidEntry(section_,Yfile))
    {
      fname  = cf_.GetStr(section_,Yfile);
      std::cerr << "Reading " << fname << std::endl;
      Y_Mapping.readVar(fname);
    }
    Y_Mapping = Y_Mapping * ARCMIN_TO_RADIANS;
    SD.addEventVar(Y_Mapping);
  }

  if(cf_.ValidEntry(section_,det_az_name))
    det_az  = cf_.GetReal(section_,det_az_name);

  if(cf_.ValidEntry(section_,det_el_name))
    det_el  = cf_.GetReal(section_,det_el_name);

  det_az = Y_Mapping.Interpol_ndx(det_az);
  det_el = X_Mapping.Interpol_ndx(det_el);
  SolarHighEl_ = X_Mapping.VInterpol_ndx(Track_Low_X_);
  SolarLowEl_  = X_Mapping.VInterpol_ndx(Track_High_X_);
  SolarHighAz_ = Y_Mapping.VInterpol_ndx(Track_High_Y_);
  SolarLowAz_  = Y_Mapping.VInterpol_ndx(Track_Low_Y_);
  for (int i=0; i<CenterOffsetsVect_.size(); i++)
    CenterOffsetsVect_[i] = X_Mapping.VInterpol_ndx(CenterOffsetsVect_[i]);
}




void SunSensorCalc::determineDetectorLockDownPosition(void)
{
  EventVar<double> SunCenterAz;

  SunCenterAz = (SolarLowAz_ + SolarHighAz_) / 2.0;
  DetectorAz_ = (SunCenterAz  - det_az) * TrueSolarExtent_ /  MeasuredExoExtentAz_;
  DetectorEl_ = (SolarHighEl_ - det_el) * TrueSolarExtent_ /  MeasuredExoExtentEl_;
}


void SunSensorCalc::determineCenterSumsElevation(void)
{

  for(int i = 0; i<CenterOffsetsVect_.size(); i++)
      CenterOffsetsVect_[i] = det_el - CenterOffsetsVect_[i];
}


void SunSensorCalc::determineSolarExtent(void)
{
  double exoExtent;
  double stddev;
  std::string str;
  EventVar<double> tmp;
  EventVar<double> timegrid;
  EventVar<double> lowel;
  SolarExtent_   = SolarHighEl_ - SolarLowEl_;
  SolarExtentAz_ = SolarHighAz_ - SolarLowAz_;
  
  valarray<bool> good  = ((zImpact_ >= 100.0) &&
			  (zImpact_ <= 140.0));
  
  tmp = SolarExtent_[good];
  MeasuredExoExtentEl_ = tmp.Mean();
  tmp = SolarExtentAz_[good];
  MeasuredExoExtentAz_ = tmp.Mean();
 
  EventVar<double> mstd;

  SolarExtent_.setName("SolarExtent");
//  SolarExtentCompare_.addEventVar(SolarExtent_);
  
  NormalizedSolarExtent_ = SolarExtent_ / MeasuredExoExtentEl_;

  if(cf_.GetFlag(section_,"PlotCompareSolarExtent"))
  {
    EventVar<double> tmpp(SolarExtent_);
    tmpp.setName("Model");
    SolarExtentCompare_.addEventVar(SolarExtent_);
    str = "Event " + GATS_Utilities::ConvertToString(eventNumber_);
//    SolarExtentCompare_.Plot("Solar Extent",str);
  }
  if(cf_.GetFlag(section_,"PlotSolarExtent"))
    SolarExtent_.Plot("SolarExtent","");
  
}


void SunSensorCalc::doDriftCorrection(void)
{
  DriftCorrection dc;

  dc.doLinearFit(SolarExtent_,zImpact_,SunSensorTimeGrid_);
  SolarExtent_ = dc.doDriftCorrection(TrueSolarExtent_);
  NormalizedSolarExtent_ = SolarExtent_ / TrueSolarExtent_;
  if(cf_.GetFlag(section_,"PlotCompareSolarExtent"))
  {
    EventVar<double> tmpp(SolarExtent_);
    tmpp.setName("DriftCorrected");
    SolarExtentCompare_.addEventVar(SolarExtent_);
    std::string str = "Event " + GATS_Utilities::ConvertToString(eventNumber_);
    SolarExtentCompare_.Plot("SolarExtent::doDriftCorrection",str);
  }
}



void SunSensorCalc::adjustToScienceTimeStamp(Event& L1)
{

  std::string str, tmpstr;
  std::stringstream oss;
  double SS_TimeShift = 0.0;
  std::string offset_name;

  if(srssflag_ == 'r')
    offset_name = "SunSensorTimeShift_SR";
  else if(srssflag_ == 's')
    offset_name = "SunSensorTimeShift_SS";
  else
    offset_name = "SunSensorTimeShift";

  if(cf_.ValidEntry(section_,offset_name))
    SS_TimeShift  = cf_.GetReal(section_,offset_name);

 
  oss << eventNumber_;
  str = "Event " + oss.str();

  EventVar<double> xloc_ts, csum_ts, hiy_ts, loy_ts, xsum_ts;


  xloc_ts = SunSensorTimeGrid_ - SS_TimeShift;
  csum_ts = SunSensorTimeGrid_ - SS_TimeShift;
  hiy_ts  = SunSensorTimeGrid_ - SS_TimeShift;
  loy_ts  = SunSensorTimeGrid_ - SS_TimeShift;
  xsum_ts = SunSensorTimeGrid_ - SS_TimeShift;

  SolarHighEl_ = SolarHighEl_.VInterpol(xloc_ts, ScienceTimeGrid_);
  SolarLowEl_  = SolarLowEl_.VInterpol (xloc_ts, ScienceTimeGrid_);
  SolarHighAz_ = SolarHighAz_.VInterpol(hiy_ts,  ScienceTimeGrid_);
  SolarLowAz_  = SolarLowAz_.VInterpol (loy_ts,  ScienceTimeGrid_);
  SolarExtent_ = SolarExtent_.VInterpol(xloc_ts, ScienceTimeGrid_);
  NormalizedSolarExtent_ = NormalizedSolarExtent_.VInterpol(xloc_ts, ScienceTimeGrid_);

  DetectorAz_  = DetectorAz_.VInterpol (xloc_ts, ScienceTimeGrid_);
  DetectorEl_  = DetectorEl_.VInterpol (xloc_ts, ScienceTimeGrid_);

  CenterAvesVect_ = CenterAvesVect_.VInterpol(csum_ts,ScienceTimeGrid_);
}


void SunSensorCalc::copyDataToEvent(Event& eventObj)
{
  initNames();
  eventObj.replaceEventVar(ScienceTimeGrid_);
  eventObj.replaceEventVar(Track_Low_X_);
  eventObj.replaceEventVar(Track_High_X_);
  eventObj.replaceEventVar(Track_Low_Y_);
  eventObj.replaceEventVar(Track_High_Y_);
  eventObj.replaceEventVar(SolarHighEl_);
  eventObj.replaceEventVar(SolarLowEl_);
  eventObj.replaceEventVar(SolarHighAz_);
  eventObj.replaceEventVar(SolarLowAz_);
  eventObj.replaceEventVar(SolarExtent_);
  eventObj.replaceEventVar(DetectorEl_);
  eventObj.replaceEventVar(DetectorAz_);
  eventObj.replaceEventVar(SunSensorTimeGrid_);
  eventObj.replaceEventVar(NormalizedSolarExtent_);
  eventObj.replaceEventVar(CenterAvesVect_);
  eventObj.replaceEventVar(CenterSumsVect_);
  eventObj.addEventVar(QualityOfFit_);
}


void SunSensorCalc::correctCenterSumsLocations()
{
  std::stringstream oss;
  std::string str, tmpstr;
  std::map<int, EventVar<double> > CenterAveMap;
  std::map<int, EventVar<double> > CenterOffsetMap;
  std::map<int, EventVar<double> >::iterator iter;
  std::vector<int> rowvect;
  int nvars = CenterAvesVect_.size();
  int len   = SunSensorTimeGrid_.size();

  valarray<bool> good  = ((SunSensorTimeGrid_ >= 70.0) &&
			  (SunSensorTimeGrid_ <= 71.0));

  EventVar<int> Rows;
  EventVar<double> Aves;
  EventVar<double> Offs;
  EventVar<double> tmp(len);

  int row;
  CenterOffsets_.resize(nvars);
  for(int j=0; j<nvars; j++)
  {
    Rows = CenterRowsVect_[j];
    Rows = Rows[good];
    row  = Rows.max();
    rowvect.push_back(row);
    str = "CenterAve_" + GATS_Utilities::ConvertToString(row);
    tmp.setName(str);
    tmp.setValueToMissing();
    CenterAveMap[row] = tmp;
    str = "CenterOffset_" + GATS_Utilities::ConvertToString(row);
    tmp.setName(str);
    tmp.setValueToMissing();
    CenterOffsetMap[row] = tmp;
  }
  for(int j=0; j<nvars; j++)
  {
    Rows = CenterRowsVect_[j];
    Aves = CenterAvesVect_[j];
    Offs = CenterOffsetsVect_[j][good];
    CenterOffsets_[j] = Offs.Mean();
    for(int i=0; i<Rows.size(); i++)
    {
      row = Rows[i];
      iter = CenterAveMap.find(row);
      if(iter != CenterAveMap.end())
      {
	CenterAveMap[row][i]    = Aves[i];
      }
    }
  }
  CenterAvesVect_.clear();
  CenterAvesVect_.setName("AveragePixelValues");
  CenterOffsets_.setName("SolarCenterSumsOffsets");
  for(int j=0; j<rowvect.size(); j++)
  {
    row  = rowvect[j];
    Aves = CenterAveMap[row];
    CenterAvesVect_.addEventVar(Aves);
  }
  
}




void SunSensorCalc::dumpData(std::string base, int grid)
{
  std::ofstream outfile;
  std::string fname = base + GATS_Utilities::ConvertToString(eventNumber_) + ".dat";
  outfile.open(fname.c_str(),std::ofstream::out);
  EventVar<double> timegrid;

  timegrid = (grid > 0) ? ScienceTimeGrid_ : SunSensorTimeGrid_;
  int len = timegrid.size();
  for(int i=0;i<len; i++)
  {
    outfile << timegrid[i] << ", ";
    outfile << Track_Low_X_[i] << ", " << Track_High_X_[i] << ", ";
    outfile << Track_Low_Y_[i] << ", " << Track_High_Y_[i] << ", ";
    outfile << Track_High_X_[i]-Track_Low_X_[i] << ", ";
    outfile << SolarLowEl_[i] << ", " << SolarHighEl_[i] << ", ";
    outfile << SolarExtent_[i] << ", " << NormalizedSolarExtent_[i] << ", ";
    outfile << DetectorEl_[i]  << ", " << DetectorAz_[i];
    for(int j = 0; j < 31; j++)
      outfile << ", " << AvePixVect_[j][i];
    outfile << "\n";
  }
  outfile.close();
}