//
//  $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:	RefrReg
//
// Filename:    RefrReg.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 "RefrReg.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;

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

RefrReg::RefrReg(Event &InEvent, Event &Stat, ConfigFile &cf, int regNum):AltReg(InEvent, Stat, cf, regNum)
//!< Constructor for RefrReg Class
//!< Loads Initial RefrReg Data from file
{
#ifdef LOGOUT
  cout << "In Constructor RefrReg::RefrReg(Event &, Event &, ConfigFile &)" << endl;
#endif
   //Initialize To Zero ... Set in setConfigParameters
   //altRange[0] = 0.0;
   //altRange[1] = 0.0;
   //gridRange[0] = 0.0;
   //gridRange[1] = 0.0;
   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);
   //z37.resize(iter, MISSINGVAL);
   setConfigParameters(cf, regNum);

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

   //allSignals.setName("allSignals");

   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");

}

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

void RefrReg::calcZImpact()
//!< Finds TP Location corresponding to Angle Impact
{
#ifdef LOGOUT
   cout << "In Function RefrReg::calcZImpact()" << endl;
#endif

   EventVar<double> z;

   z = radiusSC;
   z *= cos(angleImpact);
   z -= radiusC[0];

   putMissingValues(z, EVENTVAR_ALTMIN, EVENTVAR_ALTMAX);

   zImpact = z;

   return;

}

void RefrReg::compareAltRegMethods(EventVarVect<double> methods){
//!<  Compare Current Registration with input Registrations
#ifdef LOGOUT
cout << "In Function RefrReg::compareAltRegMethods(EventVarVec<double> methods)" << endl;
#endif

   int i;
   int evvSize = methods.size();
   double value;
   double offset;
   double band13; 
   int arrSize = methods[0].size();
   EventVar<double> method(arrSize);

   //setZcompareParams();

   cout << "CompareNumbers:  " << zCompareIndex << "  " << zCompareRatio << endl;
   //cout << "Band13 / PassThru / Band07 / Refr01 / Refr02" << endl;
   cout << "Band13 / Band07 / Refr01 / PassThru / Refr02" << endl;

   method = methods[0];
   band13 = method[zCompareIndex] + zCompareRatio * (method[zCompareIndex+1] - method[zCompareIndex]);
   offset = 40.0 - band13;

   for(i=0; i<evvSize; i++){
      method = methods[i];
      arrSize = methods[i].size();
      methodComparisonZValues[i] = method[zCompareIndex] + zCompareRatio * (method[zCompareIndex+1] - method[zCompareIndex]) + offset;
      cout << methodComparisonZValues[i] << "\t";
   }
   //cout << "Final:  "  << zImpact[zCompareIndex] << "  " << zImpact[zCompareIndex+1] << "  " << zCompareRatio << endl;
   value = zImpact[zCompareIndex] + zCompareRatio * (zImpact[zCompareIndex+1] - zImpact[zCompareIndex]) + offset;
   methodComparisonZValues[i] = value;
   cout << value << endl;
   cout << "Offset:  " << offset << endl;

   return;

}

void RefrReg::writeComparisonValues(Event &currEvent){
//!<  Compare Current Registration with input Registrations
#ifdef LOGOUT
cout << "In Function RefrReg::writeComparisonValues(Event currEvent)" << endl;
#endif
   writeToEvent(currEvent, methodComparisonZValues);
   writeToEvent(currEvent, methodComparisonBands);
   writeToEvent(currEvent, methodComparison);

  return;

}

void RefrReg::plotRefractionAngles(){
//!<  Create Refraction Angle Plots for Validation 
#ifdef LOGOUT
cout << "In Function RefrReg::plotRefractionAngles()" << endl;
#endif

   int size = angleGeomEdgeBottom.size();

   EventVar<double> angle(size);
   EventVar<double> z(size);
   EventVar<double> refr1(size);

   char fname[50];
   //EventVar<double> refr2(size);

   sprintf(fname, "SolarShrinkageOnRefrTime.out");
   solarShrinkageSolar.writeVar(fname);

   sprintf(fname, "SolarShrinkageOnDetTime.out");
   solarShrinkage.writeVar(fname);

   sprintf(fname, "RefrAngleBEdgeOnRefrTime.out");
   angleRefrSolar.writeVar(fname);

   sprintf(fname, "RefrAngleOnDetTime.out");
   angleRefr.writeVar(fname);

   angle = angleGeomEdgeBottom;
   angle -= angleRefr;

   z = radiusSC;
   z *= cos(angle);
   z -= radiusC[0];

   sprintf(fname, "GeomAngleBottom.out");
   angle.writeVar(fname);

   z.setAxisRange(-100, 300);
   z.setName("Bottom Edge Refraction");
   z.Plot("AltitudeRegistration", "Refraction Validation", angleRefr);
   //z.setName("SolarShrinkageAtBE");
   //z.Plot("AltitudeRegistration", "Refraction Validation", solarShrinkage);

   sprintf(fname, "AltitudeBEdgeOnDetTime.out");
   z.writeVar(fname);

   refr1 = angleRefr;
   refr1 -= solarShrinkage;
   angle -= normSolarExtent;
   //angle -= exoSolarExtent[0];
   z = radiusSC;
   z *= cos(angle);
   z -= radiusC[0];

   sprintf(fname, "GeomAngleTop.out");
   angle.writeVar(fname);

   sprintf(fname, "SolarExtent.out");
   normSolarExtent.writeVar(fname);

   sprintf(fname, "AltitudeTEdgeOnDetTime.out");
   z.writeVar(fname);

   sprintf(fname, "RefrAngleTEdgeOnDetTime.out");
   refr1.writeVar(fname);

   refr1.setName("Top Edge Refraction Angle");
   refr1.setDescription("Top Edge Refraction Angle");
   z.setName("Top Edge Refraction");
   z.Plot("AltitudeRegistration", "Refraction Validation", refr1);
   //z.setName("SolarShrinkageAtTE");
   //z.Plot("AltitudeRegistration", "Refraction Validation", solarShrinkage);

   //FOV Plot
   //angle += lockdownAngle;
   //z = radiusSC;
   //z *= cos(angle);
   //z -= radiusC[0];
   //z.setName("FOV Refraction");
   //z.Plot("AltitudeRegistration", "Refraction Validation - FOV", );

   return;
}

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

   int iMin, iMax;
   int holder;
   int eventSize = angleImpact.size();
   char strIt[56];
   double *impactAngleByIter;

   cout << "Refraction Registration with "  << refrRegNum << " Iterations." << endl;
   impactAngleByIter = new double[refrRegNum];

   //EventVarVec<double> zRefrCompare;
   EventVar<double> angleImpactComposite("CompositeI", "Impact Angle", "Radians", 0.0, eventSize);
   EventVar<double> zImpactComposite("CompositeZ", "Impact Altitudee", "Km", 0.0, eventSize);

   iteration = 0;

   if(logGriddedRefraction){
      cout << "Modelled Refraction" << endl;
      logGridData(S.angleRefr, S.zImpact);
   }

   if(logDebug) cout << "ExoExtentVal:  " << exoSolarExtent[0] << endl;
   interpolateRefr();

 //  checkInputs();

   normSolarExtent *= exoSolarExtent[0];

   //plotRefractionAngles();
/*
   sprintf(strIt, "Solar Extent");
   Extent.setName(strIt);
   Extent.setDescription("");
   Extent.setDescription(strIt);
   Extent.setAxisRange(0.3, 0.5);
   Extent.Plot("Refr RefrReg", strIt, altRegTime);
*/

   //angleRefr = angleRefrSolar;

   getStartIndex(S);

   if(logDebug){
      cout << "ExoAtmosphericValues: ZImpact  SolarExtent  ExoSolarExtentVal NormalizedValue" << endl;
      if(eventType == 1 && sIndex+660<normSolarExtent.size())
         cout << zImpact[sIndex+660] << "\t" << normSolarExtent[sIndex+660] << "\t" << exoSolarExtent[0] 
              << "\t" << normSolarExtent[sIndex+660]/exoSolarExtent[0] << endl;
      if(eventType == 2 && sIndex>660)
         cout << zImpact[sIndex-660] << "\t" << normSolarExtent[sIndex-660] << "\t" << exoSolarExtent[0] 
              << "\t" << normSolarExtent[sIndex-660]/exoSolarExtent[0] << endl;
   }

   holder = sIndex;

   //altRegTime.setAxisRange();
   //zImpact.setAxisRange();

   for(int i=0; i<refrRegNum; i++){

      if(logDebug) cout << "Iteration Loop sIndex:  " << sIndex << endl;

      interpolateRefr(S);

      fillAngleImpact();

      angleImpactComposite += angleImpact;

      impactAngleByIter[i] = angleImpact[holder];

      if(eventType == 1)
         sIndex--;
      if(eventType == 2)
         sIndex++;

      if(plotRefrReg){
         sprintf(strIt, "Composite It%02i", i+1);
         angleImpact.setName(strIt);
         angleImpact.setDescription("");
         angleImpact.setDescription(strIt);
         angleImpact.setAxisRange(0.3, 0.5);
         plotAnglesIt.addEventVar(angleImpact);
      }
   }

   sIndex = holder;

   angleImpactComposite /= refrRegNum;

   angleImpact = angleImpactComposite;

   calcRMS(impactAngleByIter, holder);  //Convert to Z and Output to LogFile

   delete[] impactAngleByIter;

   if(plotRefrReg){

      zGeom = radiusSC;
      zGeom *= cos(angleImpact);
      zGeom -= radiusC[0];

      sprintf(strIt, "Refraction Angle Bottom Edge");
      zGeom.setName(strIt);
      zGeom.setAxisRange(0., 100.);
      zGeom.setDescription("");
      zGeom.setDescription(strIt);
      zGeom.Plot("RefrReg:  RefraReg", strIt, angleRefr);

      sprintf(strIt, "Composite Mean");
      angleImpact.setName(strIt);
      angleImpact.setAxisRange(0.3, 0.5);
      angleImpact.setDescription("");
      angleImpact.setDescription(strIt);
      //angleImpact.Plot("Refr RefrReg", strIt, altRegTime);
      plotAnglesIt.addEventVar(angleImpact);
   }

   calcZImpact();

   angleImpactEdgeBottom = angleImpact;
   zImpactEdgeBottom = zImpact;
   zImpactEdgeBottom.setName("zImpactEdgeBottom");
   //currEvent.addEventVar(zImpactEdgeBottom);
   writeToEvent(currEvent, zImpactEdgeBottom);
   angleRefr.setName("angleRefrEdgeBottom");
   //currEvent.addEventVar(angleRefr);
   writeToEvent(currEvent, angleRefr);

   //angleImpact -= Extent[0];
   //angleImpact -= fpaSolarExtent;
   angleImpact -= normSolarExtent;

   if(plotRefrReg){

      sprintf(strIt, "Impact Angle Mean Top Edge");
      angleImpact.setName(strIt);
      angleImpact.setAxisRange(0.3, 0.5);
      angleImpact.setDescription(strIt);
      //angleImpact.Plot("Refr RefrReg", strIt, altRegTime);
      plotAnglesIt.addEventVar(angleImpact);
   }

   angleImpact += lockdownAngle;

   if(plotRefrReg){
      sprintf(strIt, "Impact Angle Mean FOV");
      angleImpact.setName(strIt);
      angleImpact.setAxisRange(0.3, 0.5);
      angleImpact.setDescription("");
      angleImpact.setDescription(strIt);
      //angleImpact.Plot("Refr RefrReg", strIt, altRegTime);
      plotAnglesIt.addEventVar(angleImpact);
      plotAnglesIt.setName("Plotting Angles");
      plotAnglesIt.setDescription("Refraction Registration Angles");
   }

   calcZImpact();
   zImpact.setDescription("RefrReg Final Impact Z FOV");

   //zImpact = radiusSC;
   //zImpact *= cos(angleImpact);
   //zImpact -= radiusC[0];

   setZcompareParams();

   plotZImpactIt.addEventVar(zImpact);
   //plotZImpactIt.setDescription("RefrReg Final Impact Z FOV");

   if(plotRefrReg){
      sprintf(strIt, "Calculated Refraction Angle");
      //zGeom.setName(strIt);
      zImpact.setAxisRange(0., 40.);
      zImpact.setDescription(strIt);
      zImpact.Plot("Altitude Registration", strIt, angleRefr);
      EventVarVect<double> Evv;
      Evv.addEventVar(fpaSolarExtent);
      Evv.addEventVar(normSolarExtent);
      Evv.Plot("Refraction Registration", "SolarExtentComparison", zImpact);
   }

   //getComparisonRange(iMin, iMax);
   iMin = sIndex;
   if(eventType == 1)
      iMax = sIndex+refrRegNum;
      //iMax = sIndex+refrRegNum+1;
   if(eventType == 2)
      iMax = sIndex-refrRegNum;
      //iMax = sIndex-refrRegNum-1;

   if(logDebug){
      cout << "iMin:  " << iMin << "iMax:  " << iMax << endl;
      cout << "CompZ: " << zImpact[iMin] << " "  << zImpact[iMax] << endl;
   }

   //Set ZOffset
   zOffset = 0.0;
   for(int i=iMin; i<=iMax; i++){
      zOffset += plotZImpactIt[1][i] - plotZImpactIt[0][i];
   }
   zOffset /= (iMax - iMin + 1);

   if(logDebug) cout << "ZOffset " << zOffset << endl;


   if(plotRefrReg){

      //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);
      //plotAnglesIt.Plot("Altitude Registration", "Lockdown Iterations", altRegTime);

      //zImpact.setAxisRange(30.0, 50.0);
      //Difference.setAxisRange(-10.0, 10.0);
      //zImpact.Plot("Altitude Registration", "Difference Plot", Difference);
      plotZImpactIt.Plot("Altitude Registration", "Refraction Registration", corrBandSignal);
      plotZImpactIt.Plot("Altitude Registration", "Refraction Registration", altRegTime);
      plotAnglesIt.Plot("Altitude Registration", "Refraction Registration", altRegTime);
   }


   cout << "Completed Refraction Registration. " << endl << endl;
//   radiusSC = radiusSC_Shifted;

   return;

}

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

   zCompareIndex = searchZImpact();  //Sets Index
   zCompareRatio = setRatio();

   if(logDebug)
      cout << "Compare Index z[ " << zCompareIndex << " ]:  " << zImpact[zCompareIndex] << "  Ratio:  " << zCompareRatio << endl;

   return;

}

void RefrReg::calcRMS(double *ang, int index)
{
#ifdef LOGOUT
cout << "In Function RefrReg::calcRMS(double*, int)" << endl;
#endif

   double sum = 0.0;
   double dif = 0.0;
   double mean = 0.0;
   double rms = 0.0;
   double *z;
   z = new double[refrRegNum];

   for(int i=0; i<refrRegNum; i++){
      //ang[i] -= fpaSolarExtent[index];
      ang[i] -= normSolarExtent[index];
      ang[i] += lockdownAngle[index];
      z[i] = radiusSC[index] * cos(ang[i]) - radiusC[0];
      sum += z[i];
   }
   mean = sum / refrRegNum;
   sum = 0.0;
   for(int i=0; i<refrRegNum; i++){
      dif = z[i] - mean;
      sum += (dif * dif);
   }
   rms = sqrt(sum / refrRegNum);

   meanSolution = mean;
   rmsSolution = rms;

   cout << "QA Refraction Registration:  "  << mean << " +/- " << rms << endl;

   delete[] z;

   return;

}

void RefrReg::getStartIndex(ModelData &M){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::getStartIndex(ModelData&)" << endl;
#endif

   int i=0;
   int foundIndex = 1;
   char sIndChar[20];
   //cout << "MaxAlt:  " << maxRefrAlt << "  AngleRefr:  " << M.angleRefr.size() << endl;
   // Scan Model for Start Altitude
   while(M.zImpact[i] > maxRefrAlt && i<M.angleRefr.size()) {
      //cout << "zImpact[ " << i << " ]:  " << M.zImpact[i] << endl;
      i++;
   }

   if(i == M.angleRefr.size()){
      cerr << "RefrReg::getStartIndex(ModelData&) Error: " << endl;
      cerr << "  Did Not Find Model Z Index: "  << maxRefrAlt;
      cerr << "                    between [ "  << M.zImpact[0] << ",  " << M.zImpact[M.zImpact.size()-1] << " ]" << endl << endl;
      throw runtime_error("Failure in AltRet::getStartIndex(ModelData&) Error: Did Not Find Initial Z Index");
   }


   //Calculate Z
   //angleImpact = angleImpactEdge;  //GeoTopEdge
   angleImpact = angleGeomEdgeBottom;  //GeoBottomEdge
   angleImpact -= angleRefr;       //ImpactBottomEdge
   angleImpact -= normSolarExtent;  //ImpactTopEdge
   angleImpact += lockdownAngle;  //ImpactFOV

   if(minLockIndex == -1 || maxLockIndex == -1)
      getLockdownSettleIndex();

   if(plotDebug){
      sprintf(sIndChar, "Impact Before");
      angleImpact.Plot("RefrReg::getStartIndex", sIndChar, altRegTime);
   }

   if(plotDebug){
      sprintf(sIndChar, "Impact After");
      angleImpact.Plot("RefrReg::getStartIndex", sIndChar, altRegTime);
   }

   calcZImpact();

   plotZImpactIt.addEventVar(zImpact);

   sIndex = angleRefr.size()/2;
   if(eventType == 1){
      while(angleRefr[sIndex] < M.angleRefr[i] && sIndex>=0) sIndex--;
      if(sIndex < 0) foundIndex = 0;
   }

   if(eventType == 2){
      while(angleRefr[sIndex] < M.angleRefr[i] && sIndex<angleRefr.size()) sIndex++;
      if(sIndex >= angleRefr.size()) foundIndex = 0;
   }

   if(foundIndex && logDebug){
      cout << "SIndex:  " << sIndex << "    Model Refr:  " << M.angleRefr[i] << "  Start Refr:  " << angleRefr[sIndex] << endl;
      cout << "SIndex:  " << sIndex << "    zImpact:  " << zImpact[sIndex] << endl;
   }

   if(plotDebug){
      sprintf(sIndChar, "SIndex %05i", sIndex);
      angleRefr.Plot("RefrReg::getStartIndex", sIndChar, altRegTime);
   }

   if(!foundIndex){
      cerr << "RefrReg::getStartIndex(ModelData&) Error: "  << endl;
      cerr << "    Did Not Find Start Refraction Index: "  << sIndex << endl;
      if(eventType == 1)
         cerr << "    EventType:  " << eventType << "  angleRefr[0]  " << angleRefr[0]
              << "  angleRefr[m]  " << angleRefr[angleRefr.size()/2]  << endl;
      else
         cerr << "    EventType:  " << eventType << "  angleRefr[n-1]  " << angleRefr[angleRefr.size()-1]
              << "  angleRefr[m]  " << angleRefr[angleRefr.size()/2]  << endl;
      throw runtime_error("Failure in AltRet::getStartIndex() Error: Did Not Find Initial Z Index");
   }

return;

}

void RefrReg::interpolateRefr(ModelData &S){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::interpolateRefr(ModelData &S)" << endl;
#endif

   double zTmp;
   int i;
   int foundRefr = 0;

   //Search Model Refraction For Indices Above and Below RefractionAngle Value at sIndex
   //Interpolate On this to get Model Impact Altitude at RefractionAngle
   //Convert Model Impact Altitude to Impact Angle at sIndex
   for(i=0; i<S.angleRefr.size()-2; i++){
      if(S.angleRefr[i] <= angleRefr[sIndex] && S.angleRefr[i+1] > angleRefr[sIndex]){
         foundRefr = i;
         interpLine(angleRefr[sIndex], S.angleRefr[i], S.angleRefr[i+1],
                    //S.angleImpact[i+1],  S.angleImpact[i], angleImpact[sIndex]);
                    S.zImpact[i],  S.zImpact[i+1], zTmp);
         angleImpact[sIndex] = acos((radiusC[0] + zTmp) / radiusSC[sIndex]);
         i=S.angleRefr.size();
      }
   }

   //Log If Found and Desired
   if(foundRefr && logDebug){
      cout << "InterpolateRefr() Found Index at " << foundRefr << " with ImpactAngle of " << angleImpact[sIndex] << endl;
      cout << "Initial Refraction Angle:  " << angleRefr[sIndex] << "between [ " << S.angleRefr[foundRefr]
           << ", " << S.angleRefr[foundRefr+1] << " ]" << endl;
      cout << "        Interpolated Z  :  " << zTmp << "between [ " << S.zImpact[foundRefr]
           << ", " << S.zImpact[foundRefr+1] << " ]" << endl;
      cout << "RadCurv:  " << radiusC[0] << " RadS/C:  " << radiusSC[sIndex] << endl << endl;
      cout << "Geometric Angle:  " << angleGeom[sIndex] << " with deltaG " << angleGeom[sIndex]-angleGeom[sIndex-1] 
           << " at Index " << sIndex << endl << endl;
      //cout << "Surrounding Geom:  " << endl;
      //for(i=sIndex-100; i<sIndex+100; i++)
      //   cout << i << " has Geometric Angle of " << angleGeom[i] << " with deltaG " << angleGeom[i]-angleGeom[i-1] << endl;
      //cout << endl;
      
   }

   //Stop RefractionRegistration if Not Found
   if(!foundRefr){
      cerr << "RefrReg::interpolateRefr(ModelData) Error: Did Not Find Initial Refraction From Model: "  << angleRefr[sIndex] << endl;
      cerr << "    Model[0]:  " << S.angleRefr[0] << "  Model[n]:  " << S.angleRefr[S.angleRefr.size()-1] << endl;
      regValidityFlag = 1;
      //throw runtime_error("Failure in RefrReg::interpolateRefr(ModelData) Error: Did Not Find Initial Refraction From Model");
   }

   return;

}

void RefrReg::interpolateRefr(){
//!< Interpolate Refraction Angle To Observed Time Grid
#ifdef LOGOUT
cout << "In Function RefrReg::interpolateRefr()" << endl;
#endif

   for(int i=0; i<altRegTime.size(); i++){
      //if(altRegTime[i] < refAngTime[0]) angleRefr[i] = angleRefrSolar[0];
      //else if(altRegTime[i] > refAngTime[refAngTime.size()-2]) angleRefr[i] = angleRefrSolar[refAngTime.size()-2];
      if(altRegTime[i] < refAngTime[0]) {
         angleRefr[i] = angleRefrSolar[0];
         solarShrinkage[i] = solarShrinkageSolar[0];
      }
      else if(altRegTime[i] > refAngTime[refAngTime.size()-2]) {
         angleRefr[i] = angleRefrSolar[refAngTime.size()-2];
         solarShrinkage[i] = solarShrinkageSolar[refAngTime.size()-2];
      }
      else {
         for(int j=0; j<refAngTime.size()-2; j++){
            if(altRegTime[i] >= refAngTime[j] && altRegTime[i] < refAngTime[j+1]){
               interpLine(altRegTime[i], refAngTime[j], refAngTime[j+1], angleRefrSolar[j], angleRefrSolar[j+1], angleRefr[i]);
               interpLine(altRegTime[i], refAngTime[j], refAngTime[j+1], solarShrinkageSolar[j], solarShrinkageSolar[j+1], solarShrinkage[i]);
               //angleRefr[i] *= 3.14159265359/(180.0*60.0);
               j = refAngTime.size();
            }
         }
      }
   }

   //angleRefr.Plot("AltitudeRegistration", "Interpolate Refraction", altRegTime);

   return;
}

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

   fillAngleImpactUp();
   fillAngleImpactDown();
   return;

}



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

   double deltaI, deltaG, deltaR;

   for(int i=sIndex; i<angleImpact.size()-2; i++){
      deltaG = angleGeomEdgeBottom[i+1] - angleGeomEdgeBottom[i];
      deltaR = angleRefr[i+1] - angleRefr[i];
      deltaI = deltaG - deltaR;
      angleImpact[i+1] = angleImpact[i] + deltaI;
      //cout << "FillUp:  " << i << "\t" << angleImpact[i] << "\t" << angleImpact[i+1]  << "\t" << angleGeomEdgeBottom[i]  << "\t" << angleRefr[i] << endl;
   }

   return;

}

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

   double deltaI, deltaG, deltaR;

   for(int i=sIndex; i>0; i--){
      deltaG = angleGeomEdgeBottom[i] - angleGeomEdgeBottom[i-1];
      deltaR = angleRefr[i] - angleRefr[i-1];
      deltaI = deltaG - deltaR;
      angleImpact[i-1] = angleImpact[i] - deltaI;
      //cout << "FillDn:  " << i << "\t" << angleImpact[i] << "\t" << angleImpact[i-1]  << "\t" << angleGeomEdgeBottom[i]  << "\t" << angleRefr[i] << endl;
   }

   return;

}

int RefrReg::searchZImpact(){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::searchZImpact()" << endl;
#endif

   int num = -1;

   if(logDebug){
      cout << "In Search Z Impact Function:  " << endl;
      cout << "  zMin[ " << minLockIndex << " ]:  " << zImpact[minLockIndex] << endl;
      cout << "  zMax[ " << maxLockIndex << " ]:  " << zImpact[maxLockIndex] << endl;
   }

   for(int i=minLockIndex; i<maxLockIndex-1; i++){
      if((zImpact[i] <= zCompareAlt && zImpact[i+1] >  zCompareAlt) ||
         (zImpact[i] >=  zCompareAlt && zImpact[i+1] < zCompareAlt)){
        //zCompareIndex = i;
        num = i;
        i = maxLockIndex;
      }
   }

   return num;

}

double RefrReg::setRatio(){
//!<  Finds the Ratio Between Indices For Predetermined Z
#ifdef LOGOUT
cout << "In Function RefrReg::setRatio()" << endl;
#endif

   double r;

   r = (zCompareAlt - zImpact[zCompareIndex]) / (zImpact[zCompareIndex+1] - zImpact[zCompareIndex]);

   return r;

}

int RefrReg::getZcompareIndex(){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::getZcompareIndex()" << endl;
#endif

   return zCompareIndex;

}

double RefrReg::getZcompareRatio(){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::getZcompareRatio()" << endl;
#endif

   return zCompareRatio;

}


void RefrReg::getRefractionQA(double& rms, double& m){
//!<
#ifdef LOGOUT
cout << "In Function RefrReg::getRefractionData(double&, double&)" << endl;
#endif

  if(altRegType != 3){
     rms = MISSINGVAL;
     m = MISSINGVAL;
     cout << "WARNING:  Trying to Obtain Refraction QA Parameter From Non-Refraction Registration Object." << endl;
  }
  else {
     rms = rmsSolution;
     m = meanSolution;
  }

   return;

}