//
//  $Id$
//-----------------------------------------------------------------------
//
//            (c) Copyright 2006 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.
//
//-----------------------------------------------------------------------
//
// Class:	CO2Reg
//
// Filename:    CO2Reg.cpp
//
// Author:      Christopher W. Brown  <c.w.brown@gats-inc.com>
//
// Date:        December 2006
//
//-----------------------------------------------------------------------
// Modification History:
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

#include "ConfigFile.h"
#include "EventVar.h"
#include "AltReg.h"
#include "CO2Reg.h"
#include "ModelData.h"
#include "NonLinearity.h"
#include "sofiefov.h"
#include "ConvFunctions.h"
#include <iostream>
#include <string>
#include <stdexcept>
#include<vector>
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//
//#define LOGOUT
#define MISSINGVAL -1.e+24
#define MINRANGE1 5.0    //Minimum Altitude Range to base Pressure Registration
#define MINRANGE2 0.02   //Minimum Transmission Range to base Pressure Registration
#define EVENTVAR_ALTMIN -200.0   //Minimum Altitude for Altitude EventVars
#define EVENTVAR_ALTMAX 500.0    //Maximum Altitude for Altitude EventVars

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

//
//-----------------------------------------------------------------------
// Utility Routines:
//-----------------------------------------------------------------------
//

using namespace std;

CO2Reg::CO2Reg():AltReg()
//!< Default Constructor for AltReg Class
{
}

CO2Reg::CO2Reg(Event &InEvent, Event &Stat, ConfigFile &cf, int regNum):AltReg(InEvent, Stat, cf, regNum)
//!< Constructor for CO2Reg Class
//!< Loads Initial CO2Reg Data from file
{
#ifdef LOGOUT
  cout << "In Constructor CO2Reg::CO2Reg(Event &, Event &, ConfigFile &)" << endl;
#endif

   //initializeParams();

   //setConfigParameters(cf, regNum);

   //initializeNames();

}

CO2Reg::~CO2Reg()
//!< Default Destructor for CO2Reg Class
{
#ifdef LOGOUT
   cout << "In Destructor CO2Reg::~CO2Reg()" << endl;
#endif
}

void CO2Reg::initializeNames()
//
{

   altRegTime.setName("Registered Time");
   corrBandSignal.setName("corrBandSignal");
   angleImpact.setName("angleImpact");
   angleImpactEdge.setName("angleImpactEdge");
   angleRefr.setName("angleRefr");
   angleGeom.setName("angleGeom");
   angleGeomEdge.setName("angleGeomEdge");
   zGeom.setName("zGeom");
   zImpact.setName("zImpact");
   zImpactEdgeTop.setName("zImpactTopEdge");
   zImpactEdgeBottom.setName("zImpactEdgeBottom");
   methodComparisonZValues.setName("altReg_ZImpact");
   methodComparisonBands.setName("altReg_Band");
   methodComparison.setName("altReg_Method");
   lockdownAngle.setName("LockdownAngle");
   //EventVar<double> losVectorSun();
   radiusSC.setName("radiusSC");
   radiusSC_Shifted.setName("radiusSC_Shifted");
   //losSun.setName("losSun");
   allSignals.setName("allSignals");

}

void CO2Reg::initializeParams()
//!< Initialize Parameters
{

   gridSpacing = 0.0;
   rangeMethod = 1;
   band = 0;
   plotGridData = 0;
   plotSimData = 0;
   initTimeOffset = 0.0;
   deltaTime = 0.0;
   dRMS = 0.0;
   convCriteria = 0.0;
   maxIterations = 0;
   logModelData = 0;
   plotLockdownFunc = 0;
   plotRefrReg = 0;
   plotFunctions = 0;
   plotInputParams = 0;
   presetGrid = 0;
   presetRange = 0;
   plotDebug = 0;
   plotGeomAngles = 0;
   plotConvolve = 0;
   plotSLDC = 0;
   plotSunSensor = 0;
   plotNonLinearity = 0;
   plotModelCalcFinal = 0;
   plotModelCalcIt = 0;
   plotDifferences = 0;
   plotDifferencesIt = 0;
   plotCorrectionBand = 0;
   plotZImpact = 0;
   doConvolution = 0;
   doNonLinearCorr = 0;
   doSmoothLockdown = 0;
   logDebug = 0;
   logGriddedRefraction = 0;
   //Initialize To Defaults
   currConverge = 100.0; //Set High so it will run
   saveConverge = 100.0; //Set High so it will run
   eventType = 0;
   gridNum = 0;
   iteration = 0;
   numBands = 24;
   sIndex = 0;
   zCompareAlt = 0.0;
   zCompareRatio = 0.0;
   zCompareIndex = 0;
   timeOffset = 0.0;
   saveTimeOffset = 0.0;
   minTimeOffset = -3.0;
   //minTransmission = 0.3;
   //maxTransmission = 0.4;
   maxTimeOffset = 3.0;
   minLockIndex = -1;
   maxLockIndex = -1;
   MAXFINALTIMEOFFSET = 2.0;
   regValidityFlag = 0;
   altRange.resize(2, 0.0);
   gridRange.resize(2, 0.0);

   return;

}

void CO2Reg::calcImpactData(ModelData& S)
//!< Adds Interpolated Refraction Angle to angleGeomEdge
//!< From Model Refraction Relationship Between thetaImpact and zImpact
//!< Input(s): CO2Reg populated with Sim Data (Later CO2RegModel)
{

   int modelSize = S.angleRefr.size();
   double modelMax, modelMin, modelMaxRefr;
   double arMax, arMin, arNew;
   double aiMax, aiMin, aiNew;
   double agMax, agMin;
   EventVar<double> newAngleR(angleGeomEdge.size());
   EventVar<double> newAngleI(angleGeomEdge.size());
   EventVar<double> newAngleIM(angleRefr.size());
   EventVar<double> newAngleIB(angleGeomEdge.size());
   EventVar<double> newAngleRB(angleRefr.size());

   for(int i=0; i<angleGeomEdge.size(); i++){
      //SetUpModel
      S.angleImpact = acos((radiusC[0] + S.zImpact) / radiusSC[i]);
      S.angleGeom = S.angleImpact;
      S.angleGeom += S.angleRefr;
      if(plotModelCalcIt && i==600){
         EventVarVect<double> Evv;
         char title[35];
         sprintf(title, "calcImpactData(CO2Reg) LOOP: %04i", i);
         S.angleImpact.setAxisRange();
         S.angleRefr.setAxisRange();
         S.angleGeom.setAxisRange();
         Evv.addEventVar(S.angleImpact);
         Evv.addEventVar(S.angleRefr);
         Evv.addEventVar(S.angleGeom);
         S.zImpact.setAxisRange();
         Evv.Plot("Altitude Registration",title,S.zImpact);
      }

      modelMax = S.angleGeom[modelSize-1];
      modelMin = S.angleGeom[0];
      modelMaxRefr = S.angleRefr[modelSize-1];

      //TopEdgeCalculation
       if(angleGeomEdge[i] < modelMin){
         newAngleI[i] = angleGeomEdge[i];
         newAngleR[i] = 0.0;
       }
       else if(angleGeomEdge[i] > modelMax){
          newAngleI[i] = angleGeomEdge[i] - modelMaxRefr;
          newAngleR[i] = modelMaxRefr;
       }
       else {
          for(int j=0; j<modelSize-1; j++){
             agMin = S.angleGeom[j];
             agMax = S.angleGeom[j+1];
             if(angleGeomEdge[i] < agMax && angleGeomEdge[i] >= agMin){
                
                //Interpolate On Impact Angle
                aiMin = S.angleImpact[j];
                aiMax = S.angleImpact[j+1];
                interpLine(angleGeomEdge[i], agMin, agMax, aiMin, aiMax, aiNew);
                //Interpolate On Refraction Angle
                //arMin = S.angleRefr[j];
                //arMax = S.angleRefr[j+1];
                //interpLine(angleGeom[i], agMin, agMax, arMin, arMax, arNew);
                //cout << "Ang[" << S.zImpact[j] << "] (G, R, I, G-R):  " << angleGeom[i] << ", " << arNew << ", " << aiNew << ", " << angleGeom[i] - arNew << endl;
                newAngleI[i] = aiNew;
                //newAngleR[i] = arNew;
                newAngleR[i] = angleGeomEdge[i] - aiNew;
                j = modelSize+1;
             }
          }
       }

   if(angleGeomEdgeBottom[i] < modelMin){
      newAngleIB[i] = angleGeomEdgeBottom[i];
      newAngleRB[i] = 0.0;
   }
   else if(angleGeomEdgeBottom[i] > modelMax){
      newAngleIB[i] = angleGeomEdgeBottom[i] - modelMaxRefr;
      newAngleRB[i] = modelMaxRefr;
   }
   else {
      for(int j=0; j<modelSize-1; j++){
         agMin = S.angleGeom[j];
         agMax = S.angleGeom[j+1];
         if(angleGeomEdgeBottom[i] < agMax && angleGeomEdgeBottom[i] >= agMin){
             //Interpolate On Impact Angle
             aiMin = S.angleImpact[j];
             aiMax = S.angleImpact[j+1];
             interpLine(angleGeomEdgeBottom[i], agMin, agMax, aiMin, aiMax, aiNew);
             //Interpolate On Refraction Angle
             //arMin = S.angleRefr[j];
             //arMax = S.angleRefr[j+1];
             //interpLine(angleGeom[i], agMin, agMax, arMin, arMax, arNew);
             //cout << "Ang[" << S.zImpact[j] << "] (G, R, I, G-R):  " << angleGeom[i] << ", " << arNew << ", " << aiNew << ", " << angleGeom[i] - arNew << endl;
             newAngleIB[i] = aiNew;
             //newAngleR[i] = arNew;
             newAngleRB[i] = angleGeomEdgeBottom[i] - aiNew;
             j = modelSize+1;
          }
       }
    }
}

   angleImpactEdge = newAngleI;
   angleRefr = newAngleR;
   angleImpactEdgeBottom = newAngleIB;
   angleRefrEdgeBottom = newAngleRB;

   //angleGeomEdge = angleImpactEdge + angleRefr;

   if(plotModelCalcFinal){
      EventVarVect<double> eVV;
      char title[30];
      sprintf(title,"calcImpactData(CO2Reg)");
      //evv.addEventVar(angleImpact);
      eVV.addEventVar(angleImpactEdge);
      eVV.addEventVar(angleRefr);
      eVV.addEventVar(angleGeomEdge);
      eVV.addEventVar(angleImpactEdgeBottom);
      eVV.addEventVar(angleRefrEdgeBottom);
      eVV.addEventVar(angleGeomEdgeBottom);

      altRegTime.setAxisRange();
      eVV.Plot("Altitude Registration",title,altRegTime);

   }

   return;

}

void CO2Reg::readAttenSettings(Event& InEvent)
//!< Reads Attenuator Setting EventVar from Event Variable
//!< Attenuator Setting Located In Level 0 Event Object
{

   InEvent.getEventVar("L0_AttenuatorSettings", attenSetting);

   return;

}

/*
void CO2Reg::readGeometricData(Event& InEvent, Event& SD, ConfigFile& CF)
//!< Reads Geometric EventVars from Event Variable
//!< Reads Config File for Registration Parameters
//!< Stores Information from TimeReg and SigCorr in CO2Reg Class Object
//!< Input(s):
//!<      Event with EventVars SimSig, SimAlt, SimRef
//!< Stored Output(s): 
{
#ifdef LOGOUT
   cout << "In Function CO2Reg::readGeometricData(Event& InEvent, Event& SD, ConfigFile& CF)" << endl;
#endif

   EventVar<double> obsSize;
   InEvent.getEventVar("Final_Signal_Time", obsSize);

//Resize ALL EventVars
   corrBandSignal.resize(obsSize.size());
   altRegTime.resize(obsSize.size());
   angleImpact.resize(obsSize.size());
   angleImpactEdge.resize(obsSize.size());
   angleRefr.resize(obsSize.size());
   angleGeom.resize(obsSize.size());
   angleGeomEdge.resize(obsSize.size());
   zGeom.resize(obsSize.size());
   zImpact.resize(obsSize.size());
   lockdownAngle.resize(obsSize.size());
   radiusSC.resize(obsSize.size());
   radiusSC_Shifted.resize(obsSize.size());
   attenSetting.resize(numBands);
   signalExo.resize(numBands);
   //losSun.resize(obsSize.size());

   angleRefrSolar.resize(obsSize.size());
   lockdownAngleAz.resize(obsSize.size());
   fpaSolarExtent.resize(obsSize.size());
   normSolarExtent.resize(obsSize.size());
   exoSolarExtent.resize(1);

   //Get EventVars From InEvent
   cout << "Band:  " << band+1 << endl;
   InEvent.getEventVar("Corrected_Signals", allSignals);
   corrBandSignal = allSignals[band];
   corrBandSignal.setName("Altitude_Registration_Band");

   InEvent.getEventVar("Final_Signal_Time", altRegTime);
   altRegTime.setName("L1_AltRegTime");

   InEvent.getEventVar("L1_TPAlt", zGeom);

   zGeom.setName("Geometric Center Altitude Z");
   zGeom.setAxisRange(0.0, 180.0);

   InEvent.getEventVar("DetectorLockdown_EL", lockdownAngle);
   lockdownAngle.setName("L1_AngleLockdownEl");

   InEvent.getEventVar("DetectorLockdown_AZ", lockdownAngleAz);
   lockdownAngleAz.setName("L1_AngleLockdownAz");

   InEvent.getEventVar("SolarHighElevationEdge", solarHighEdgeEl);
   solarHighEdgeEl.setName("SolarHighElevationEdge");

   InEvent.getEventVar("SolarLowElevationEdge", solarLowEdgeEl);
   solarLowEdgeEl.setName("SolarLowElevationEdge");

   InEvent.getEventVar("SolarHighAzimuthEdge", solarHighEdgeAz);
   solarHighEdgeAz.setName("SolarHighAzimuthEdge");

   InEvent.getEventVar("SolarLowAzimuthEdge", solarLowEdgeAz);
   solarLowEdgeAz.setName("SolarLowAzimuthEdge");

   InEvent.getEventVar("SolarExtent", fpaSolarExtent);
   fpaSolarExtent.setName("SolarExtent");

   InEvent.getEventVar("L1_SolarExtent", exoSolarExtent);
   exoSolarExtent.setName("L1_SolarExtent");

   InEvent.getEventVar("NormalizedSolarExtent", normSolarExtent);
   normSolarExtent.setName("NormalizedSolarExtent");

   InEvent.getEventVar("L1_SCRad", radiusSC);
   radiusSC.setName("L1_SCRad");

   InEvent.getEventVar("L1_DopplerShift", doppler);
   doppler.setName("L1_DopplerV");

   InEvent.getEventVar("L1_CurvatureRadius", radiusC);
   radiusC.setName("L1_CurvatureRadius");

   InEvent.getEventVar("L1_Drift_Signal_Ref", signalExo);
   signalExo.setName("L1_VExo");

   InEvent.getEventVar("L1_RefractionAngles", angleRefrSolar);
   //InEvent.getEventVar("L1_RefractionAngles", angleRefrSolar);
   //angleRefrSolar.setName("L1_VExo");
   //angleRefrSolar *= 3.14159265359/180.0;

   //InEvent.getEventVar("SolarTrackerTimes", refAngTime);
   InEvent.getEventVar("L1_RefractionTime", refAngTime);
   //refAngTime.setName("L1_SolarTrackTime");
   refAngTime.setName("L1_RefractionTime");

   //InEvent.getEventVar("L1_Drift_Time_Ref", exoTime);
   //exoTime.setName("L1_timeVExo");

   //InEvent.getEventVar("L1_SolarLineOfSight", losSun);
   //losSun.setName("L1_SolarLineOfSight");

   //Initialize Remaining EventVars
   angleImpact.setValueConst(MISSINGVAL);
   angleImpact.setName("L1_AngleImpact");
   angleImpact.setUnits("Radians");
   angleImpact.setDescription("Calculated Modeled FOV Impact Angle");

   angleImpactEdge.setValueConst(MISSINGVAL);
   angleImpactEdge.setName("L1_AngleImpactEdge");
   angleImpactEdge.setUnits("Radians");
   angleImpactEdge.setDescription("Calculated Observed Edge Impact Angle");

   angleRefr.setValueConst(MISSINGVAL);
   angleRefr.setName("L1_RefractionAngles");
   angleRefr.setUnits("Radians");
   angleRefr.setDescription("Calculated Modeled Refraction Angle");

   radiusSC_Shifted.setValueConst(MISSINGVAL);
   radiusSC_Shifted.setName("L1_SCRad_TimeShifted");

   angleGeom.setValueConst(MISSINGVAL);
   angleGeom.setName("L1_AngleGeom");

   angleGeomEdge.setValueConst(MISSINGVAL);
   angleGeomEdge.setName("Tmp_AngleGeomEdge");

   zImpact.setValueConst(MISSINGVAL);
   zImpact.setName("L1_SimulatedAngleGeom");

   //Store Event Type
   if(zGeom[2] > zGeom[1]) eventType = 1;
   else eventType = 2;

   //signalExo = 1;
   //signalExo = signalExoVar[band];
   //for(int i=0; i<corrBandSignal.size()-1; i++)
   //  if(zGeom[i] < 153.0 && zGeom[i] > 152.0){
   //     signalExo = corrBandSignal[i];
   //     i = corrBandSignal.size()+1;
   //}
   corrBandSignal = corrBandSignal / signalExo[band];
   cout << "EventType:  " << eventType << endl;
   if(logDebug) cout << "Exo Var:  " << signalExo[band] << endl;

   if(plotInputParams && band == 12){
      string title;
      //title = "Interpolate Angle " + iteration;
      altRegTime.setAxisRange();
      altRegTime.setDescription("Observation Time");
      altRegTime.setUnits("Sec");

      corrBandSignal.setAxisRange();
      corrBandSignal.setDescription("Observation Transmission");
      corrBandSignal.setUnits("Trans");
      title = "Observed Band Input";
      corrBandSignal.Plot("Altitude Registration Input",title,altRegTime);

      zGeom.setAxisRange();
      title = "Observed Z Input";
      zGeom.Plot("Altitude Registration Input",title,altRegTime);

      zGeom.setAxisRange();
      title = "Observed Correction Band";
      corrBandSignal.Plot("Altitude Registration Input",title,zGeom);

      angleRefrSolar.setAxisRange();
      refAngTime.setAxisRange();
      title = "Refraction Angle";
      angleRefrSolar.Plot("Altitude Registration Input",title,refAngTime);
      //angleRefrSolar.Plot("Altitude Registration Input",title,altRegTime);

      lockdownAngle.setAxisRange();
      title = "Observed Lockdown Input";
      lockdownAngle.Plot("Altitude Registration Input",title,altRegTime);

      lockdownAngleAz.setAxisRange();
      title = "Observed Lockdown Azimuthal Input";
      lockdownAngleAz.Plot("Altitude Registration Input",title,altRegTime);

      fpaSolarExtent.setAxisRange();
      title = "FPA Solar Extent Input";
      fpaSolarExtent.Plot("Altitude Registration Input",title,altRegTime);

      normSolarExtent.setAxisRange();
      title = "FPA Normalized Solar Extent Input";
      normSolarExtent.Plot("Altitude Registration Input",title,altRegTime);

      solarHighEdgeEl.setAxisRange();
      title = "FPA Solar High Elevation Edge Input";
      solarHighEdgeEl.Plot("Altitude Registration Input",title,altRegTime);

      solarLowEdgeEl.setAxisRange();
      title = "FPA Solar Low Elevation Edge Input";
      solarLowEdgeEl.Plot("Altitude Registration Input",title,altRegTime);

      solarHighEdgeAz.setAxisRange();
      title = "FPA Solar High Azimuthal Edge Input";
      solarHighEdgeAz.Plot("Altitude Registration Input",title,altRegTime);

      solarLowEdgeAz.setAxisRange();
      title = "FPA Solar Low Azimuthal Edge Input";
      solarLowEdgeAz.Plot("Altitude Registration Input",title,altRegTime);

      radiusSC.setAxisRange();
      title = "Observed S/C Radius Input";
      radiusSC.Plot("Altitude Registration Input",title,altRegTime);

      doppler.setAxisRange();
      title = "Observed Doppler Input";
      doppler.Plot("Altitude Registration Input",title,altRegTime);

      title = "Observed Correction Band Input";
      zGeom.Plot("Altitude Registration Input",title,corrBandSignal);

      title = "Observed Lockdown Angle Input";
      zGeom.Plot("Altitude Registration Input",title,lockdownAngle);

      //title = "Observed/Measured Refraction Angle Input";
      //zGeom.Plot("Altitude Registration Input",title,angleRefr);

      title = "Observed SCRadius Input";
      zGeom.Plot("Altitude Registration Input",title,radiusSC);

      title = "Observed Doppler Input";
      zGeom.Plot("Altitude Registration Input",title,doppler);

   }

}

void CO2Reg::setConfigParameters(ConfigFile& C, int regNum)
//!< Reads Config File for Registration Parameters
//!< Stores Information into CO2Reg Class Object Private Variables
{
#ifdef LOGOUT
   cout << "In Function CO2Reg::setConfigParameters(ConfigFile& C)" << endl;
#endif

    vector<int> temp;


   //Otherwise uses default Values used in Constructor
   try {
      //band = C.GetInt("AltitudeRegistration", "RegistrationBand");
      temp  = C.GetIntValues("AltitudeRegistration", "RegistrationBand");
      band = temp[regNum];
      band--;
      b37.resize(temp.size());
      for(int i = 0; i<temp.size(); i++)
         b37[i] = double (temp[i]);
         //b37[i] = temp[i];
   }
   catch(...)
   {
      band = 12;
      b37.resize(1);
      b37[0] = double (band);
      //b37[0] = band;
      cout << "Registration Band Set to Default Value:  " << band << endl;
   }

   try {
      temp  = C.GetIntValues("AltitudeRegistration", "AltRegType");
      //temp  = C.GetRealValues("AltitudeRegistration", "AltRegType");
      altRegType = temp[regNum];
      //altRegType = C.GetInt("AltitudeRegistration", "AltRegType");
      m37.resize(temp.size());
      for(int i = 0; i<temp.size(); i++)
         m37[i] = double (temp[i]);
         //m37[i] = temp[i];

   }
   catch(...)
   {
      altRegType = 2;
      m37.resize(1);
      m37[0] = double (altRegType);
      cout << "AltRegType Set to Default:  " << altRegType << endl;
   }

   try {
      plotGridData = C.GetInt("AltitudeRegistration", "PlotGriddedEventVars");
   }
   catch(...)
   {
      plotGridData = 0;
   }

   try {
      plotDebug = C.GetInt("AltitudeRegistration", "PlotDebug");
   }
   catch(...)
   {
      plotDebug = 0;
   }

   try {
      fovFileName = C.GetStr("AltitudeRegistration", "FOVFile");
   }
   catch(...)
   {
      sprintf(fovFileName, "../Data/FOV/sofie_fov_V1.1_elevation.txt"); 
   }

   try {
      plotGeomAngles = C.GetInt("AltitudeRegistration", "PlotGeometricAngles");
   }
   catch(...)
   {
      plotGeomAngles = 0;
   }

   try {
      plotConvolve = C.GetInt("AltitudeRegistration", "PlotConvolution");
   }
   catch(...)
   {
      plotConvolve = 0;
   }

   try {
      plotSLDC = C.GetInt("AltitudeRegistration", "PlotSLDC");
   }
   catch(...)
   {
      plotSLDC = 0;
   }

   try {
      plotModelCalcIt = C.GetInt("AltitudeRegistration", "PlotIteratedModelCalc");
   }
   catch(...)
   {
      plotModelCalcIt = 0;
   }

   try {
      plotModelCalcFinal = C.GetInt("AltitudeRegistration", "PlotModelCalc");
   }
   catch(...)
   {
      plotModelCalcFinal = 0;
   }

   try {
      plotDifferences = C.GetInt("AltitudeRegistration", "PlotDifferences");
   }
   catch(...)
   {
      plotDifferences = 0;
   }

   try {
      plotDifferencesIt = C.GetInt("AltitudeRegistration", "PlotIteratedDifferences");
   }
   catch(...)
   {
      plotDifferencesIt = 0;
   }

   try {
      plotCorrectionBand = C.GetInt("AltitudeRegistration", "PlotBandSignal");
   }
   catch(...)
   {
      plotCorrectionBand = 0;
   }

   try {
      plotAllBands = C.GetInt("AltitudeRegistration", "PlotAllBands");
   }
   catch(...)
   {
      plotAllBands = 0;
   }

   try {
      plotZImpact = C.GetInt("AltitudeRegistration", "PlotZ");
   }
   catch(...)
   {
      plotZImpact = 0;
   }

   try {
      doConvolution = C.GetInt("AltitudeRegistration", "ConvolveFOV");
   }
   catch(...)
   {
      doConvolution = 1;
   }

   try {
      doNonLinearCorr = C.GetInt("AltitudeRegistration", "AddNonLinearity");
   }
   catch(...)
   {
      doNonLinearCorr = 1;
   }

   try {
      doSmoothLockdown = C.GetInt("AltitudeRegistration", "SmoothLockdown");
   }
   catch(...)
   {
      doSmoothLockdown = 0;
   }

   try {
      logDebug = C.GetInt("AltitudeRegistration", "LogDebug");
   }
   catch(...)
   {
      logDebug = 0;
   }

   try {
      logGriddedRefraction = C.GetInt("AltitudeRegistration", "LogRefraction");
   }
   catch(...)
   {
      logGriddedRefraction = 0;
   }

   try {
      refrRegNum = C.GetInt("AltitudeRegistration", "RefractionProfiles");
   }
   catch(...)
   {
      refrRegNum = 5;
   }

   try {
      plotSimData = C.GetInt("AltitudeRegistration", "PlotModel");
   }
   catch(...)
   {
      plotSimData = 0;
   }

   try {
      initTimeOffset = C.GetReal("AltitudeRegistration", "Init_Time_Offset");
   }
   catch(...)
   {
      initTimeOffset = 0.02;
      cout << "Initial TimeOffset Set to Default Value:  " << initTimeOffset << endl;
   }

   try {
      maxRefrAlt = C.GetReal("AltitudeRegistration", "RefractionMaxAlt");
   }
   catch(...)
   {
      maxRefrAlt = 30.0;
      cout << "Maximum Refraction Altitude Set to Default Value:  " << maxRefrAlt << endl;
   }

   try {
      zCompareAlt = C.GetReal("AltitudeRegistration", "ComparisonAltitude");
   }
   catch(...)
   {
      zCompareAlt = 50.0;
      cout << "Comparison Altitude Set to Default Value:  " << zCompareAlt << endl;
   }

   try {
      MAXFINALTIMEOFFSET = C.GetReal("AltitudeRegistration", "Maximum_Time_Offset");
   }
   catch(...)
   {
      MAXFINALTIMEOFFSET = 2.0;
      cout << "Maximum Final TimeOffset Set to Default Value:  " << MAXFINALTIMEOFFSET << endl;
   }

   try {
      convCriteria = C.GetReal("AltitudeRegistration", "Convergence_Criteria");
   }
   catch(...)
   {
      convCriteria = 0.1;
      cout << "Convergence Criteria Set to Default Value:  " << convCriteria << endl;
   }
      dRMS = convCriteria + 5.0;

   try {
      maxIterations = C.GetInt("AltitudeRegistration", "Maximum_Iterations");
   }
   catch(...)
   {
      maxIterations = 10;
      cout << "Maximum Iterations Set to Default Value:  " <<  maxIterations << endl;
   }

   try {
      logModelData = C.GetInt("AltitudeRegistration", "LogCreateModelGrid");
   }
   catch(...)
   {
      logModelData = 0;
   }

   try {
      plotFunctions = C.GetInt("AltitudeRegistration", "PlotFunctionIO");
   }
   catch(...)
   {
      plotFunctions = 0;
   }

   try {
      plotLockdownFunc = C.GetInt("AltitudeRegistration", "PlotLockdownFunction");
   }
   catch(...)
   {
      plotLockdownFunc = 0;
   }

   try {
      plotRefrReg = C.GetInt("AltitudeRegistration", "PlotRefractionRegistration");
   }
   catch(...)
   {
      plotRefrReg = 0;
   }

   try {
      plotSunSensor = C.GetInt("AltitudeRegistration", "PlotSunSensor");
   }
   catch(...)
   {
      plotSunSensor = 0;
   }

   try {
      plotNonLinearity = C.GetInt("AltitudeRegistration", "PlotNonLinearity");
   }
   catch(...)
   {
      plotNonLinearity = 0;
   }

   try {
      plotInputParams = C.GetInt("AltitudeRegistration", "PlotInputParameters");
   }
   catch(...)
   {
      plotInputParams = 0;
   }

   try {
      gridSpacing = C.GetReal("AltitudeRegistration", "GridSpacingInKM");
   }
   catch(...)
   {
      gridSpacing = 0.2;
      cout << "Grid Spacing Set to Default Value:  " <<  gridSpacing << endl;
   }
//Added 20070803
   try {
      gridRange[0] = C.GetReal("AltitudeRegistration", "GridRangeLower");
      gridRange[1] = C.GetReal("AltitudeRegistration", "GridRangeUpper");
      if(gridRange[0] >= gridRange[1] - gridSpacing) throw;
      presetGrid = 1;
   }
   catch(...)
   {
      gridRange[0] = 3.0;
      gridRange[1] = 180.0;
      presetGrid = 0;
      cout << "Grid Range Set to Default Value:  [ " << gridRange[0] << ", " << gridRange[1]  << endl;
   }

   try {
      //rangeMethod = C.GetInt("AltitudeRegistration", "ComparisonMethod");
      altRange[0] = C.GetReal("AltitudeRegistration", "ComparisonRangeLower");
      altRange[1] = C.GetReal("AltitudeRegistration", "ComparisonRangeUpper");
      if(altRange[0] < 2.0 && altRange[1] < 2.0) rangeMethod = 2;
      else rangeMethod = 1;

      if(rangeMethod == 1 && (altRange[0] >= altRange[1] - MINRANGE1)){
         cout << "WARNING:  Improper CO2Reg Altitude Range:  { " << altRange[0] << ", " << altRange[1] << " }" << endl;
         throw;
      }
      if(rangeMethod == 2 && (altRange[0] >= altRange[1] - MINRANGE2)){
         cout << "WARNING:  Improper CO2Reg Transmission Range:  { " << altRange[0] << ", " << altRange[1] << " }" << endl;
         cout << "          Will Use Default Altitude Range." << endl;
         throw;
      }
      presetRange = 1;
   }
   catch(...)
   {
      rangeMethod = 1;
      altRange[0] = 35.0;
      altRange[1] = 45.0;
      presetRange = 0;
      cout << "Altitude Comparison Range Set to Default Value:  [ " << altRange[0] << ", " << altRange[1]  << " ]" << endl;
   }

   if(band >= numBands){
      //band = 12;
      cout << "Illegal Registration Band: " << band+1 << ".  There are only " << numBands << " Bands.  Right?" << endl;
      band = 12;
   }

   if(fabs(initTimeOffset) >= 1.0){
      cout << "That's a Large Time Offset: " << initTimeOffset << ".  ";
      initTimeOffset = 0.5;
      cout << "Resetting to " <<  initTimeOffset << " seconds.  Please keep Initial Offset Less than 1.0 seconds." << endl;
      //initTimeOffset = 0.0;
   }

   if(gridSpacing != 0)
     gridNum = int (((gridRange[1] - gridRange[0]) / gridSpacing ) + 1.01);
   else gridNum = 0;

   //Need to Figure Out how to store Range in Config File and Retrieve
   //Until Then Use default Values used in SetConfig
   //altRange = C.getRealArray("Altitude Registration", "Registration Range", num); //Might have to assign later
   //gridRange = C.getRealArray("Altitude Registration", "Level1 Grid Range", num); //Might have to assign later
   //gridSpacing = C.getRealArray("Altitude Registration", "Level1 Grid Spacing", num); //Might have to assign later
   //int num = 2;

   return;

}
*/

void CO2Reg::shiftSCRad(EventVar<double> &SC){

   int shift = getShift(2);
   double arscEdge, rscMin, rscMax;
   double tTarget, tMin, tMax;

   //Get Shifted Radius
   for(int i=minLockIndex; i<maxLockIndex; i++){
   //for(int i=0; i<altRegTime.size()-1; i++){
      //If timeOffset is OOR, then do linear extrapolate
      //Check to see not within Observed Window
      if(altRegTime[i] + timeOffset < altRegTime[minLockIndex])
         SC[i] = radiusSC[i] + timeOffset * (radiusSC[i+1] - radiusSC[i]) / (altRegTime[i+1] - altRegTime[i]);
      else if(altRegTime[i] + timeOffset > altRegTime[altRegTime.size()-1])
         SC[i] = radiusSC[i] - timeOffset * (radiusSC[i-1] - radiusSC[i]) / (altRegTime[i-1] - altRegTime[i]);
         //cout << "AngTimeOOR:  " << (altRegTime[i]+timeOffset)*1000 << "  altTime0 " << altRegTime[0] << "  altTimeN " << altRegTime[altRegTime.size()-1] << "  timeOffset " << timeOffset << endl;
      else { // Interpolate and Add Lockdown
         //shift = getShift(i);  // Find Indices to Shift
         tTarget = altRegTime[i] + timeOffset;
         tMin =  altRegTime[i+shift];
         tMax = altRegTime[i+shift+1];
         rscMin = radiusSC[i+shift];
         rscMax = radiusSC[i+shift+1];
         interpLine(tTarget, tMin, tMax, rscMin, rscMax, arscEdge);
         SC[i] = arscEdge;
         //cout << "AngTime:  " << (altRegTime[i]+timeOffset)*1000 << "  angleGeom " << angleGeom[i+shift] << "  angleImpact " << angleImpact[i] << "  angleGeom1 " << angleGeom[i+shift+1] << endl;
      }
   }
   //Assign Linear adjustment to Last radius Shifted Value
   SC[altRegTime.size()-1] = SC[altRegTime.size()-2]*2.0 - SC[altRegTime.size()-3] ;
    return;
}

void CO2Reg::calcZImpactFOV()
//!< Finds TP Location corresponding to AngleImpactFOV
{
#ifdef LOGOUT
   cout << "In Function CO2Reg::calcZImpactFOV()" << endl;
#endif

   //int i;
   int shift = getShift(2);
   double arscEdge, rscMin, rscMax;
   double tTarget, tMin, tMax;
   EventVar<double> radSCshift(radiusSC.size());

   shiftSCRad(radSCshift);

   zImpact = radSCshift;
   zImpact *= cos(angleImpact);
   zImpact -= radiusC[0];
   radiusSC_Shifted = radSCshift;

   // Verifying by Greg
   cout << "Verification of impact input relationships from calcZimpactFOV() " << ((radiusC[0] + zImpact[700])/cos(angleImpact[700])) << " should = " << radiusSC_Shifted[700] << endl;

//   cout << "MinLockIndex:  " << minLockIndex << "  Z:  " << zImpact[minLockIndex] << endl;
//   cout << "MaxLockIndex:  " << maxLockIndex << "  Z:  " << zImpact[maxLockIndex] << endl;

//   if(plotFunctions){
/*
      string title("Function:  calcZImpactFOV(Rc)");
      zImpact.setAxisRange();
      angleImpact.setAxisRange(15.0, 25.0);
      angleImpact.Plot("Altitude Registration",title,zImpact);

      zImpact.setAxisRange();
      altRegTime.setAxisRange();
      zImpact.Plot("Altitude Registration",title,altRegTime);

      zImpact.setAxisRange();
      radiusSC_Shifted.setAxisRange(6925.0, 7000.0);
      zImpact.Plot("Altitude Registration",title,radiusSC_Shifted);

      EventVarVect<double> Comparison;
      Comparison.addEventVar(radiusSC_Shifted);
      Comparison.addEventVar(radiusSC);
      Comparison.Plot("Altitude Registration", title, zImpact);
      Comparison.Plot("Altitude Registration", title, altRegTime);
*/
      char strIt[56];
      sprintf(strIt, "ImpactAltitudeZ_FOV_%02i", iteration);
      zImpact.setName(strIt);
      sprintf(strIt, "Impact Altitude Z FOV Iteration %02i", iteration);
      zImpact.setDescription(strIt);
      //nameEV = "Model Signal Iteration " + strIt;
      //signalMod.setName(nameEV);
      zImpact.setAxisRange(0.0, 180.0);

      plotZImpactIt.addEventVar(zImpact);

//   }

   return;

}

double CO2Reg::getTimeOffset()
//!< Returns value stored in Altreg.timeOffset
{
   return timeOffset;
}

void CO2Reg::convertSignals2Model(ModelData& S)
//!< Interpolate Values from Observed onto Model Grid
{
#ifdef LOGOUT
   cout << "In Function CO2Reg::convertSignals2Model(ModelData& S)" << endl;
#endif
   double pctg;
   int obsSize = corrBandSignal.size();
   int modSize = S.corrBandSignal.size();

   EventVar<double> zImpactMod(obsSize);
   EventVar<double> signalMod(obsSize);
   double oSig, osigMin, osigMax;
   double mZTarget, oZMin, oZMax;

   for(int i=0; i<obsSize; i++){
      mZTarget = zImpact[i];
      zImpactMod[i] = mZTarget;
      if(mZTarget > S.zImpact[0])  //Z Outside of Range ... Fill with Highest Model Value
         signalMod[i] = S.corrBandSignal[0];
      else if(mZTarget < S.zImpact[modSize-1])  //Z Outside of Range ... Fill with Lowest Model Value
         signalMod[i] = S.corrBandSignal[modSize-1];
      else {
         for(int j=0; j<modSize-1; j++){
            oZMin = S.zImpact[j+1];
            oZMax = S.zImpact[j];
            osigMin = S.corrBandSignal[j+1];
            osigMax = S.corrBandSignal[j];
            if(mZTarget >= oZMin && mZTarget < oZMax){
               interpLine(mZTarget, oZMin, oZMax, osigMin, osigMax, oSig);
               signalMod[i] = oSig;
      //         zImpactMod[i] = mZTarget;
               //Ey[i] = Ey[i] / Ey[0];
               //cout << "SSVals:  " << Cx[j] << "  " << Cx[j+1] << " " << Ex[i] << "  " << Cy[j] << "  " << Cy[j+1] << " " << Ey[i] << endl;
               j = modSize + 2;
            }
         }
      }
   }

   //currZImpact = zImpactObs;
   modelSignal = signalMod;

   if(plotFunctions || plotCorrectionBand){

      char strIt[56];
      sprintf(strIt, "Initial Model Signal Iteration %02i", iteration);
      signalMod.setName(strIt);
      signalMod.setAxisRange(0.0, 1.0);
      plotCBSignalIt.addEventVar(signalMod);

   }

   return;
}

void CO2Reg::initializeCO2Reg(){
//!< Initialize Values 
#ifdef LOGOUT
cout << "In Function CO2Reg::initializeCO2Reg()" << endl;
#endif
 
   cout << "CO2 Registration using Band " << band+1 << endl;

   //Initialize Values
   saveTimeOffset = timeOffset;
   iteration = 0;
   angleImpact = angleImpactEdge;
   angleGeom = angleGeomEdge;

   //zGeom.setDescription("Geometric Edge");
   plotZImpactIt.addEventVar(zGeom);

   if(plotFunctions){
      angleImpactEdge.setAxisRange(0.3, 0.5);
      plotLockdownIt.addEventVar(angleImpactEdge);
      plotLockdownIt.setDescription("Lockdown Iteration Plot");
      corrBandSignal.setAxisRange(0.1, 0.5);
      plotCBSignalIt.addEventVar(corrBandSignal);
      plotCBSignalIt.setDescription("Correction Band Transmissions");
   }
   else {
      if(plotLockdownFunc){
         angleImpactEdge.setAxisRange(0.3, 0.5);
         plotLockdownIt.addEventVar(angleImpactEdge);
         plotLockdownIt.setDescription("Lockdown Iteration Plot");
      }
      if(plotCorrectionBand){
         corrBandSignal.setAxisRange(0.1, 0.5);
         plotCBSignalIt.addEventVar(corrBandSignal);
         plotCBSignalIt.setDescription("Correction Band Transmissions");
      }
   }

   getLockdownSettleIndex();

   return;

}

void CO2Reg::doPressureRegistration(Event &currEvent, ModelData& S){
//!< Loop through Time to Match Geo with Sim
//!< Sets Time Offset For Geo Data
//!<
#ifdef LOGOUT
cout << "In Function CO2Reg::doPressureRegistration(Event& currEvent, ModelData& Sim)" << endl;
#endif

   int iMin = 0;
   int iMax = 0;

   initializeCO2Reg();

   while(fabs(dRMS) > convCriteria && iteration < maxIterations){

      if(logDebug) cout << "Newton Method Iteration:  "  << iteration << endl;

      angleImpact = angleImpactEdge;
      angleGeom = angleGeomEdge;

      addLockdown();

      calcZImpactFOV();

      if(!iteration) writePassThruToEvent(currEvent);

      convertSignals2Model(S);

      if(doConvolution)
         convolveFOV(currEvent);

      if(doNonLinearCorr)
         addNonlinearity();

      compareSignals(iMin, iMax);
      //if(iteration < 3){
      //   getComparisonRange(iMin, iMax);
      //}
      //calcRMS(iMin, iMax);

      setTimeOffset();

      //Plot Iterated Differences
      if(plotDifferencesIt)
         plotDiff();

      iteration++;
   }

   finalizeCO2Reg(iMin, iMax);

   writeMetrics(currEvent);

   return;

}


void CO2Reg::writeMetrics(Event& currEvent)
//!< Writes the Output of All Member Variables to Event Object Tmp or Level1
//!< Requires No Parameters ... Uses Config for Output Parameters and Locations
{
#ifdef LOGOUT
   cout << "In Function AltReg::writeMetrics(Event& currEvent)" << endl;
#endif

   char name[128];
   char description[128];
   cout << band << endl;
   sprintf(name, "AltReg_ShiftedTime_Band%02i", band); //PressReg_ShiftedTime
   EventVar<double> tOffsetPR(name, "Shifted Time of Pressure Registration", timeOffset, 1);
   sprintf(name, "AltReg_ShiftedAltZ_Band%02i", band); //PressReg_ShiftedZ
   EventVar<double> zOffsetPR(name, "Shifted Z of Pressure Registration", zOffset, 1);
   sprintf(name, "AltReg_RMSofFit_Band%02i", band); //PressReg_RMS
   EventVar<double> rmsPR(name, "RMS of Pressure Registration", rmsSolution, 1);
   sprintf(name, "AltReg_Method_Band%02i", band); //PressReg_Method
   EventVar<double> regMethod("PressReg_Method", "Registration Method", altRegType, 1);

   cout << "Final Time Offset is " << tOffsetPR[0] << endl;
   cout << "Final RMS is         " << rmsPR[0] << endl;

   writeToEvent(currEvent, tOffsetPR);   
   writeToEvent(currEvent, zOffsetPR);   
   writeToEvent(currEvent, rmsPR);   
   writeToEvent(currEvent, regMethod);   

/*
      name = tOffsetPR.getName();
      descript = tOffsetPR.getDescription();
      cout << "Adding " << name << ":  " << descript << endl;
      currEvent.addEventVar(tOffsetPR);

      name = zOffsetPR.getName();
      descript = zOffsetPR.getDescription();
      cout << "Adding " << name << ":  " << descript << endl;
      currEvent.addEventVar(zOffsetPR);

      name = rmsPR.getName();
      descript = rmsPR.getDescription();
      cout << "Adding " << name << ":  " << descript << endl;
      currEvent.addEventVar(rmsPR);
*/

return;
}




void CO2Reg::finalizeCO2Reg(int ia, int iz)
{
#ifdef LOGOUT
cout << "In Function CO2Reg::finalizeCO2Reg(int, int)" << endl;
#endif
//!< Writes Log Message, Plots Data, and Sets Final Values


   cout << "Altitude Registration Converged to:  " << timeOffset << " seconds at RMS:  "
        << currConverge << " in " << iteration-1 << " iterations" << endl;
   cout << "Final Altitude Comparison Range: ( " << zImpact[ia] << ", " << zImpact[iz] << " )" << endl << endl;
   //Set ZOffset
   if(iteration > 1)
      setZOffset(ia, iz);
   else cout << "ERROR:  CO2 Registration Finalized At 1 Iteration." << endl;

   if(plotFunctions){

      altRegTime.setAxisRange(20.0, 80.0);
      plotLockdownIt.Plot("Altitude Registration", "Lockdown Iterations", altRegTime);
      plotCBSignalIt.Plot("Altitude Registration", "Model Signal Iterations", altRegTime);
      zImpact.setAxisRange(0.0, 180.0);
      plotCBSignalIt.Plot("Altitude Registration", "Model Signal Iterations", zImpact);
      plotZImpactIt.Plot("Altitude Registration", "Impact Altitude Iterations", altRegTime);
      plotDiff();
   }
   else {
      if(plotLockdownFunc){
         altRegTime.setAxisRange(20.0, 80.0);
         plotLockdownIt.Plot("Altitude Registration", "Lockdown Iterations", altRegTime);
      }
      if(plotCorrectionBand){
         int plotIt = iteration*2;
         altRegTime.setAxisRange(20.0, 80.0);
         plotCBSignalIt.Plot("Altitude Registration", "Model Signal Iterations", altRegTime);
         zImpact.setAxisRange(0.0, 180.0);
         plotCBSignalIt.Plot("Altitude Registration", "Model Signal Iterations", zImpact);
         zImpact.Plot("Altitude Registration", "Model Signal Iterations", plotCBSignalIt[plotIt]);
      }
      if(plotZImpact){
         altRegTime.setAxisRange(20.0, 80.0);
         plotZImpactIt.Plot("Altitude Registration", "Impact Altitude Iterations", altRegTime);
      }
      if(plotDifferences){
         plotDiff();
      }
   }

   radiusSC = radiusSC_Shifted;

   return;

}

/*
void CO2Reg::plotDiff(){
#ifdef LOGOUT
cout << "In Function CO2Reg::plotDiff()" << endl;
#endif

   EventVar<double> Difference(altRegTime.size());
   char diffSetName[256];
   char subtitle[256];
   sprintf(diffSetName, "Iteration %02i -- ( %6.3f secs )", iteration, timeOffset);
   sprintf(subtitle, "Difference Band %02i (Model - Observed)", band+1);
   Difference.setName(diffSetName);
   Difference = modelSignal;
   Difference -= corrBandSignal;
   Difference /= modelSignal;
   Difference *= 100.0;
   zImpact.setAxisRange(30.0, 50.0);
   Difference.setAxisRange(-10.0, 10.0);
   zImpact.Plot("Altitude Registration", subtitle, Difference);
   sprintf(subtitle, "CO2 Registration Model %02i with Non-Linearity", band+1);
   zImpact.Plot("Altitude Registration", subtitle, modelSignal);

   return;
}
*/

void CO2Reg::setZOffset(int a, int z)
{
#ifdef LOGOUT
cout << "In Function CO2Reg::setZOffset(int, int)" << endl;
#endif

   zOffset = 0.0;
   for(int i=a; i<=z; i++){
      zOffset += plotZImpactIt[iteration-1][i] - plotZImpactIt[1][i];
   }
   zOffset /= (z - a + 1);

   return;

}

/*void AltReg::checkInputs(){

   getLockdownSettleIndex();

   checkSolarExtent();
   
   checkRefractionAngles();

   return;

}

void AltReg::checkSolarExtent(){

   EventVar<double> Ev(normSolarExtent);
   int lo, hi;
   int i;

   searchValidIndices(Ev, lo, hi);  // Finds Missing Value Indexes
   
   //i = lo;
   while((Ev[lo] < 0.0 || Ev[lo] > 1.0) && lo < hi) lo++;
   i = lo;
   while(Ev[i] >=0.0 && Ev[i] <= 1.0 && i < hi) i++;
   if(i < hi) hi = i;


   if(minLockIndex < lo){
      cout << "Solar Extent Indices Changed MinIndex From " << minLockIndex << " to " << lo << endl; 
      minLockIndex = lo;
   }
   if(maxLockIndex > hi){
      cout << "Solar Extent Indices Changed MaxIndex From " << maxLockIndex << " to " << hi << endl; 
      maxLockIndex = hi;
   }
   
   return;

}

void AltReg::checkRefractionAngles()
{
   EventVar<double> Ev(angleRefr);

   if(eventType == 1)


*/

void CO2Reg::convolveFOV(Event &currEvent)
//!<  Convolves FOV and SLDC with Model Data
//!<  Read FOV Data From File
//!<  Reads SLDC From EventVar
//!<  Runs Fortran Convolution Function
{
#ifdef LOGOUT
cout << "In Function CO2Reg::convolveFOV()" << endl;
#endif

   double simExo, lockdownExo;
   double bandOffset;
   char title[50];

   EventVarVect<double> Evv;
   EventVar<double> Ev;
   ifstream ifs(fovFileName, ios::in);
   //vector<double> simExo;
   vector<double> respFOV;
   //vector<double> angleFOVarcMin;
   vector<double> angleFOV;
   vector<double> respSLDC;
   vector<double> angleSLDC;
   vector<double> respModel;
   vector<double> angleImpactModel;
   vector<double> angleGeomModel;
   vector<double> angleLock;
   vector<double> respConvolved;

   int numFOVs, numObs, numSLDCs;
   int iBeg, iEnd;

   if(eventType == 1)
      bandOffset = bandOffset_SR;
   else
      bandOffset = bandOffset_SS;

   cout << "Using bandOffset Value:  "  << bandOffset << "  for band " << band+1 << endl;

   sofiefov gFOVdata;

   //string fovFileName("../Data/FOV/sofie_fov_V1.1_elevation.txt");

   iBeg = minLockIndex;
   iEnd = maxLockIndex;
   //Read FOV Data
   //sofiefov fovSOFIE;
   //Check that the fovFile exists
   if(ifs)
      gFOVdata.read_fov_file(fovFileName, 1);   /// add reverse=1 so it will reverse FOV file
   else
   {
      cerr << "  Unable To Read FOV File " << fovFileName << endl << endl;
      throw runtime_error("Could Not Read FOV File in convolveFOV().");
   }
   
   respFOV = gFOVdata.get_fov_response(band+1);
   //angleFOV = gFOVdata.get_fov_elev(band+1);
   angleFOV = gFOVdata.get_fov_elev2(band+1, bandOffset);
   numFOVs = angleFOV.size();

   //Input File is In ArcMins --- Change to Radians
   for(int i=0; i<numFOVs; i++)
      angleFOV[i] /= (60.0 * 180.0 /  3.14159265358979);

   //Read SLDC Data
   currEvent.getEventVar("L1_SolarModel", Evv);
   //Ev = Evv[band+1];
   Ev = Evv[band];
   respSLDC = vector<double>(&Ev[0], &Ev[0]+Ev.size() ); 
   currEvent.getEventVar("L1_SLDCAngleFromTopEdge", Ev);
   angleSLDC = vector<double>(&Ev[0], &Ev[0]+Ev.size() );
   numSLDCs = angleSLDC.size();
 
   if(plotSLDC && iteration == 2){
      EventVar<double> EvSLDC;
      EventVar<double> EvSLDCBase;
      EvSLDC = Evv[band];
      EvSLDCBase = Evv[2];

      cout << "Plot SLDC" << endl;
      EvSLDC.Plot("AltReg", "SLDC Band", Ev);
      EvSLDCBase.Plot("AltReg", "SLDC Baseline", Ev);
   }

/*   for(int i=0; i<numFOVs; i++)
      cout << i << " angleFOV:  " << angleFOV[i] << "  fovSignal:  " << respFOV[i] << endl;
   for(int i=0; i<numSLDCs; i++)
      cout << i << " angleSLDC:  " << angleSLDC[i] << "  sldcSignal:  " << respSLDC[i] << endl;
*/
   //Prepare Model Data
   respModel = vector<double>(&modelSignal[0], &modelSignal[0]+modelSignal.size());
   angleImpactModel = vector<double>(&angleImpact[0], &angleImpact[0]+angleImpact.size());
   angleGeomModel = vector<double>(&angleGeom[0], &angleGeom[0]+angleGeom.size());
   angleLock = vector<double>(&lockdownAngle[0], &lockdownAngle[0]+lockdownAngle.size());
   respConvolved.resize(respModel.size());
   numObs = respModel.size();

   lockdownExo = getExoLockdown();
 
   if(eventType == 1){
      invertVector(respModel);
      invertVector(angleImpactModel);
      invertVector(angleGeomModel);
      invertVector(angleLock);
   }

   setConvolutionIndices(angleImpactModel, angleGeomModel, iBeg, iEnd);
 
   ::calc_sim_exo__(&respFOV[0], &angleFOV[0], &numFOVs, &respSLDC[0], &angleSLDC[0], &numSLDCs, &lockdownExo, &simExo);
   if(iteration == 2){
      cout << "  Lockdown Exo:  " << lockdownExo << endl;
      cout << "  Response Size: " << respConvolved.size() << endl;
      cout << "  Num Observs  :  " << numObs << endl;
      cout << "  Num FOV      :  " << numFOVs << endl;
      cout << "  Num SLDC     :  " << numSLDCs << endl;
      cout << "  Sim Exo     :  " << simExo << endl;
   }
//   for(int i=0; i<respConvolved.size(); i++)
 //     cout << i << " angleGeom:  " << angleGeomModel[i] << "  angleLock:  " << angleLock[i] << " modelSignal:  " << respModel[i] << "  Angle:  "  << angleImpactModel[i] << " angleRefr:  " << angleGeomModel[i] - angleImpactModel[i] << endl;

   ::conv_fov_atm_sun__(&numObs, &iBeg, &iEnd, &angleImpactModel[0], &angleGeomModel[0], &angleLock[0], &respModel[0],
                        &respFOV[0], &angleFOV[0], &numFOVs, &respSLDC[0], &angleSLDC[0], &numSLDCs, &simExo, &respConvolved[0]); 

   normalize(respConvolved);  //Does not actually normalize .... just removes NANs

//   for(int i=iBeg; i<iEnd; i++)
//      cout << i << " modelSignal:  " << respConvolved[i] << "  Angle:  "  << angleImpactModel[i] << endl;

   if(eventType == 1){
      invertVector(respConvolved);
      invertVector(respModel);
      //invertVector(angleImpactModel);
      //invertVector(angleGeomModel);
      //invertVector(angleLock);
   }

   EventVar<double> Ev2("L1_ConvolvedModeledSignal", "Model Transmission Convolved With FOV & SLDC", &respConvolved[0], respConvolved.size());
   modelSignal = Ev2;
 
   //Ev2.setAxisRange();
   modelSignal.setAxisRange(0.0, 1.0);
   
   if(plotConvolve){  
 
      EventVar<double> Ev1("L1_OriginalModeledSignal", "Original Model Transmission", &respModel[0], respModel.size());
      EventVarVect<double> Evv("ConvolutionModelPlotVector", "Vector To Plot Signals Before and After Convolution");
      double initVal = 1.0;
      double max0, max1;
      EventVar<double> modelPlot(initVal, int(modelSignal.size()));
      EventVar<double> modelPlot0(initVal, int(Ev1.size()));
      max0 = Ev1.max();
      max1 = modelSignal.max();
      cout << "SimExo, MaxIn, MaxOut" << endl;
      cout << simExo << ", " << max0 << ", " << max1 << endl;
      modelPlot -= modelSignal;
      modelPlot0 -= Ev1;
      Evv.addEventVar(modelPlot0);
      Evv.addEventVar(modelPlot);
      sprintf(title, "Model Signal After Convolution %02i", iteration);
      zImpact.Plot("AltReg", title, modelPlot);
      sprintf(title, "Model Signal Convolution Before/After %02i", iteration);
      Evv.Plot("AltReg", title, zImpact);

   }

   return;

}

void CO2Reg::normalize(vector<double> &resp)
//!<  Remove Nan and Do Not Normalize By Exo -- Change Name 
{
#ifdef LOGOUT
cout << "In Function CO2Reg::normalize(vector<double>&)" << endl;
#endif

   int s, e, i;
   vector<double> zImpactVec;
   vector<int> badIndices;
   vector<double> badValues;
   double max = -1.0;

   zImpactVec = vector<double>(&zImpact[0], &zImpact[0]+zImpact.size() );
   if(eventType == 1)
      invertVector(zImpactVec);

    
   for(int i=0; i<resp.size(); i++)
      if(resp[i] < 0.0 || resp[i] > 2.0){
         badIndices.push_back(i);
         badValues.push_back(resp[i]);
         resp[i] = MISSINGVAL;
      }
    
   if(badIndices.size() > 0){
      s=resp.size()-1;
      e=badIndices.size()-1;
      while(badIndices[e] == s && e>0) {s--; e--;}
      s=0;
      while(badIndices[s] == s && s<badIndices.size()-1) s++;
      if(s<=e) cout << "WARNING - BAD INDICES IN CONVOLVED MODEL ";
      for(int i=s; i<=e; i++)
         cout << ", " << badIndices[i];
      cout << endl;
      for(int i=s; i<=e; i++)
         cout << ", " << badValues[i];
      cout << endl << endl;
   }
//   for(int i=0; i<resp.size(); i++)
//      if(resp[i] >  max && zImpactVec[i] < 150.0 && zImpactVec[i] > 140.0)
//         max = resp[i];

//   for(int i=0; i<resp.size(); i++){
//      if(resp[i] == max) cout << "Normalize Convolved Vector By " << max << " at " << zImpactVec[i] << " Km." << endl;
//      if(resp[i] > 0.0)
//         resp[i] /= max;
//   }

   return;

}



void CO2Reg::setConvolutionIndices(vector<double> &app, vector<double> &geo, int &s, int &e)
//!<  Impact Angle Must be Monotonic and Increasing.
{
#ifdef LOGOUT
cout << "In Function CO2Reg::setConvolutionIndices(vector<double>&, vector<double>&, int&, int&)" << endl;
#endif

   vector<double> zImpactVec;
   int minL, maxL;
   minL = minLockIndex;
   maxL = maxLockIndex;

   zImpactVec = vector<double>(&zImpact[0], &zImpact[0]+zImpact.size() );
   if(eventType == 1){ 
      invertVector(zImpactVec);
      minL = app.size() - maxLockIndex - 1;
      maxL = app.size() - minLockIndex - 1;
   }

   int i=minL;
   while(app[i] == MISSINGVAL && i < maxL) i++;
   if(logDebug) cout << "Convolution Indices:  ";
   while(geo[i] - app[i] < 0.0001 && i < maxL){
      //cout << i << "  Geo-App  " << geo[i] - app[i] << "  Geo:  " << geo[i] << "  App:  " << app[i] << endl;
      i++;
   }
   if(logDebug) cout << i << "Z " << zImpactVec[i];
   while(app[i] > app[i-1] && app[i-1] != MISSINGVAL && i>minL) i--;
   s = ++i;
   if(logDebug)  cout << ";  " << i << "Z " << zImpactVec[i] << endl;
   while(app[i] < app[i+1] && app[i+1] != MISSINGVAL && i<maxL) i++;
   e = i-15;
   
   //cout << "Convolution Indices:  " << zImpactVec[s] << ", " <<  zImpactVec[e] << endl;

   return;

}
   

void CO2Reg::invertVector(vector<double>& iV)
{
#ifdef LOGOUT
cout << "In Function CO2Reg::invertVector(vector<double>&)" << endl;
#endif

   int vSize = iV.size();
   vector<double> V(vSize);
   
   V = iV;
   for(int i=0; i<vSize; i++)
      iV[i] = V[vSize-i-1];

   return;

}


void CO2Reg::addNonlinearity(){
//!<
#ifdef LOGOUT
cout << "In Function CO2Reg::addNonlinearity()" << endl;
#endif

   int bandNum = band + 1;
   double gainNL = attenSetting[band];
   double sExo = signalExo[band];
   double trans;
   EventVar<double> nlTrans(modelSignal.size());

   nlTrans = modelSignal;

   for(int i=minLockIndex; i<maxLockIndex; i++){
      trans = modelSignal[i];
      if(trans >= 0.0 && trans <= 1.0)
         modelSignal[i] = sofie_sim_nonlin_tran(bandNum, trans, gainNL, sExo);
      else modelSignal[i] = MISSINGVAL;
   }

   if(plotNonLinearity){
      char setString[20];

      EventVar<double> nlTransPost(modelSignal.size());
      EventVar<double> nlDiff(modelSignal.size());
      nlTransPost = modelSignal;
      nlDiff = nlTransPost;
      nlDiff -= nlTrans;
      nlDiff /= nlTransPost;
      nlDiff *= 100.0;

      EventVarVect<double> nonLinCompare;
      nlTrans.setName("Uncorrected Transmission");
      nonLinCompare.addEventVar(nlTrans);
      nlTransPost.setName("NL Corrected Transmission");
      nonLinCompare.addEventVar(nlTransPost);
      sprintf(setString, "NL Trans Plot# %02i", iteration);
      nlDiff.setName(setString);
      nonLinCompare.Plot("Altitude Registration", setString, zImpact);

      sprintf(setString, "NL Diff It# %02i", iteration);
      nlDiff.setName(setString);
      sprintf(setString, "NL Diff Plot# %02i", iteration);
      nlDiff.setAxisRange();
      zImpact.Plot("Altitude Registration", setString, nlDiff);

   }

   if(plotFunctions || plotCorrectionBand){
      char strIt[56];
      sprintf(strIt, "Model Signal Iteration %02i", iteration);
      modelSignal.setName(strIt);
      modelSignal.setAxisRange(0.0, 1.0);
      plotCBSignalIt.addEventVar(modelSignal);
   }


   return;

}


void CO2Reg::compareSignals(int &iMin, int &iMax){
//!<
#ifdef LOGOUT
cout << "In Function CO2Reg::compareSignals(int &, int &)" << endl;
#endif

   cout << "Verification of impact input relationships from compareSignals() " << ((radiusC[0] + zImpact[700])/cos(angleImpact[700])) << " should = " << radiusSC_Shifted[700] << endl;

   if(iteration < 3)
      getComparisonRange(iMin, iMax);

   calcRMS(iMin, iMax);

   return;

}

double CO2Reg::getZPTcompare(int idx, double rat){
//!<  Returns Pass Through Altitude From CO2 Registration
#ifdef LOGOUT
cout << "In Function CO2Reg::getZPTcompare(index, ratio)" << endl;
#endif

   double ptc;

   if(altRegType != 2){
      cout << "AltRegType:  " << altRegType << endl;
      throw runtime_error("Failure in AltReg::getZPTcompare(int):  Attempting to Retrieve Uncalculated Z Value.");
   }

   ptc = plotZImpactIt[1][idx] + rat * (plotZImpactIt[1][idx+1] - plotZImpactIt[1][idx]);

   return ptc;

}

void CO2Reg::getComparisonRange(int& min, int& max){
//!< Find Start/End Indices corresponding to Altitude Range or Transmission Range for Comparison
#ifdef LOGOUT
cout << "In Function CO2Reg::getComparisonRange()" << endl;
#endif

   int i;

   EventVar<double> ComparisonVar(zImpact.size());

   if(rangeMethod == 2 && (altRange[0] > 1.0 || altRange[1] > 1.0)){
      cerr << "CO2Reg::getComparisonRange(int& min, int& max) Error: range does not match method" << endl;
      cerr << "Method:  " << rangeMethod << " Range:  ( " << altRange[0] << ", " << altRange[1] << " )" << endl;
      throw runtime_error("Failure in CO2Reg::getComparisonRange(int& min, int& max) Error: range does not match method");
   }
   if(rangeMethod != 2 && (altRange[0] < 1.1 || altRange[1] < 1.1)){
      cerr << "CO2Reg::getComparisonRange(int& min, int& max) Error: range does not match method" << endl;
      cerr << "Method:  " << rangeMethod << " Range:  ( " << altRange[0] << ", " << altRange[1] << " )" << endl;
      throw runtime_error("Failure in CO2Reg::getComparisonRange(int& min, int& max) Error: range does not match method");
   }

   if(rangeMethod == 2) ComparisonVar = corrBandSignal;
   else ComparisonVar = zImpact;

   if(minLockIndex == -1)  //Check to See if Lockdown Settle has been run
      //i = zImpact.size()/2;
      getLockdownSettleIndex();

   i=minLockIndex;

   min = -1;
   max = -1;

   if(eventType == 1){
      while(ComparisonVar[i] < altRange[0]) i++;
      min = i;
      while(ComparisonVar[i] < altRange[1]) i++;
      max = i-1;
      //cout << "Comparison Index[ " << min << ", " << max << " ]:  { " << zImpact[min] << " , " << zImpact[max] << " }; " << endl;
   }
   else if(eventType == 2){
      while(ComparisonVar[i] > altRange[1]) i++;
      min = i;
      while(ComparisonVar[i] > altRange[0]) i++;
      max = i-1;
      //cout << "Comparison Index[ " << min << ", " << max << " ]:  { " << zImpact[min] << " , " << zImpact[max] << " }; " << endl;
   }
   else return;

   return;

}

void CO2Reg::setTimeOffset(){
//!< Calculate New Time Offset If Not Converged
//!< Sets Time Offset
//!<
#ifdef LOGOUT
cout << "In Function CO2Reg::setTimeOffset()" << endl;
#endif

   // Calculate Raphson-Newton Xn+1   RMS    deltaRMS / delta t0 -- Enter Print Statement Here For Evaluation
   if(iteration == 0){   //Set default value
      timeOffset = initTimeOffset;
      deltaTime = initTimeOffset;
      saveConverge = currConverge;
   }
   //else if(fabs(currConverge) > convCriteria && iteration < maxIterations){
   else if(fabs(dRMS) > convCriteria && iteration < maxIterations){
      if(logDebug){
         cout << "LastConverge: " << saveConverge << " DRMS: " << currConverge - saveConverge << "  DTime:  "  << deltaTime << endl;
         cout << "MinTimeOffset:  " << minTimeOffset << " MaxTimeOffset: " << maxTimeOffset << endl;
      }
      timeOffset = timeOffset - currConverge / ((currConverge - saveConverge) / deltaTime);
      if(logDebug)  cout << "CalculatedTimeOffset:  " << timeOffset;
      if(timeOffset == saveTimeOffset || fabs(timeOffset - saveTimeOffset) < 1e-18){
         cout << "Newton Raphson provided Same Output Time Offset without Convergence.  Something Amiss?  " << endl;
         iteration = maxIterations + 7;
      }
      else if(timeOffset < -30.0 || timeOffset > 30.0){
         cout << "Newton Raphson provided Large Time OFfset.  Divergence?" << endl;
         iteration = maxIterations + 7;
      }
      //else if(((currConverge - saveConverge) / (timeOffset - saveTimeOffset)) > 9999999.9
      else {
         //deltaTime = timeOffset - saveTimeOffset;
         dRMS =  fabs(currConverge - saveConverge);
         if(timeOffset > saveTimeOffset) minTimeOffset = saveTimeOffset;
         if(timeOffset < saveTimeOffset) maxTimeOffset = saveTimeOffset;
         if(timeOffset < minTimeOffset) timeOffset = (maxTimeOffset + minTimeOffset) / 2.0;
         if(timeOffset > maxTimeOffset) timeOffset = (maxTimeOffset + minTimeOffset) / 2.0;
         if(logDebug) cout << "   UsingTimeOffset:  " << timeOffset << endl <<endl;
         if(fabs(timeOffset) > MAXFINALTIMEOFFSET)
            cout << "WARNING:  TIMEOFFSET is Out Of Range."  << endl;
         deltaTime = timeOffset - saveTimeOffset;
         saveConverge = currConverge;
         saveTimeOffset = timeOffset;
      }
   }
   else
      cout << "Altitude Registration Converged to:  " << timeOffset << " seconds at RMS:  " << currConverge << endl;

   return;

}

void CO2Reg::calcRMS(int iLow, int iHi)
//!<
{
   int num = iHi - iLow + 1;
   double sum = 0.0;
   double sum2 = 0.0;

   for(int i=iLow; i<=iHi; i++){
      //cout << "RMSValues:  " << sig1[i]  << "  " <<  sig2[i] << endl;
      sum += corrBandSignal[i] - modelSignal[i];                                             //Using dSignal
//      sum += (corrBandSignal[i] - modelSignal[i]) * (corrBandSignal[i] - modelSignal[i]);  //Using dRMS
//      sum2 += (corrBandSignal[i] - modelSignal[i]);                                        //Using dRMS
      sum2 += (corrBandSignal[i] - modelSignal[i]) * (corrBandSignal[i] - modelSignal[i]);   //Using dSignal and outputting RMS

      //else
      //cout << "RMSValues:  " << z1[i] << "  " <<  altRange[0]  << "  " <<  z2[i] << "  " <<  altRange[1] << endl;
   }
//   if(sum2 == 0.0) sum2 = 1.0; //sum2 is only used for the direction of convergence ...  so set to positive for 0
//   currConverge = sqrt(sum/num) * fabs(sum2)/sum2;
   currConverge = sum;
   rmsSolution = sqrt(sum2/num);
   //cout << "TimeOffset:  "  << timeOffset << "  Convergence:  sqrt(" << sum << " / " <<  num  << ") = " <<  currConverge << endl;

   return;
}