/***************************************************************************
 *   (C) Copyright 2006 by GATS, Inc.
 *         11864 Canon Blvd., STE 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: Satellite
 *
 * Filename: Satellite.cpp
 *
 * Description: Satellite has member functions to manipulate date events
 *              and process SOFIE Level1 Time Registration for a single event
 *
 * Author:      C.W.Brown
 *              (757)952-1048
 *              (757)873-5920
 *              <c.w.brown@gats-inc.com>
 *
 ***************************************************************************/
 
 #include <iostream>
 #include <fstream>
 //#include <math.h>
 #include <cmath>
 #include "Ephemeris.h"
 #include "Satellite.h"
 #include "GeoCal.h"
 #include "DVector.h"
 #include "EventVar.h"
 #include "GeneralFunctions.h"

 #define DEBUG
 //#define MAKEPLOTS 
 //#define OUTPUTLLA

 using namespace std;

 const double SUNRAD = 695500.0;


Satellite::Satellite():Ephemeris()
//!< Default Constructor for Satellite Class
{
}

Satellite::Satellite(ifstream &file):Ephemeris(file)
//!< Constructor for Satellite Class 
//!< Loads Initial Ephemeris from file
{}

Satellite::Satellite(int size)
//Satellite::Satellite(int size):Ephemeris(size){}
//!< Constructs an empty SatelliteClass object
//!< Allocates memory for 'size' data points in the stateVector
//!< This constructor will be used with the Satellite Ephemeris
{
  numPts = size;
  stateVector = new StateVector[size];
  elevation = new double[size];          //!< Elevation Angle Array for Each Time Step  (in degrees)
  tpLat = new double[size];              //!< Geodetic Tangent Point Latitude for Each Time Step (in degrees)
  tpLon = new double[size];              //!< Geodetic Tangent Point Longitude for Each Time Step (in degrees)
  tpAlt = new double[size];              //!< Geodetic Tangent Point Altitude for Each Time Step (in km)
  scLat = new double[size];              //!< Geodetic Spacecraft Latitude for Each Time Step (in degrees)
  scLon = new double[size];              //!< Geodetic Spacecraft Longitude for Each Time Step (in degrees)
  scAlt = new double[size];              //!< Geodetic Spacecraft Altitude for Each Time Step (in km)
  doppler = new double[size];            //!< Doppler Shift for Each Time Step
  sinkRate = new double[size];           //!< Elevation Angle Velocity for Each Time Step
  objectDistance = new double[size];     //!< Distance From Satellite to Object (in Km)
  earthSunDistance = new double[size];     //!< Distance From Satellite to Object (in Km)
  scRad = new double[size];          //!< Distance From Satellite to Object (in Km)
  //radCurvature = new double[size];     //!< Instantenous Radius of Curvature (in km)
  //losVector = new LOSVector[size];     //!< Instantenous Radius of Curvature (in km)
  losVectors = new DVector*[size];
  losSunVectors = new DVector*[size];
  for(int i=0; i<size; i++) {
     losVectors[i] = new DVector(3);   
     losSunVectors[i] = new DVector(3);   
  } 
  solarRadius = -1.0e+24;

//!< Instantenous Radius of Curvature (in km)
//for(int i=0; i<size; i++)
//	losVector[i].los = new DVector(3);
//tangentXYZ[6];

}

Satellite::~Satellite()
//!< Default Destructor for Satellite Class
{
}

//Satellite::Satellite:Ephemeris(mysql Database)
//!< Constructor for Satellite Class 
//!< Loads Initial Ephemeris from Database
//!< EventVar etime must be in Modified Julian Day 
//!< MJD = Julian Day - 2400000.5 

void Satellite::getEphemeris(Ephemeris& ephem, EventVar <double> etime)
//!< Function initializes Satellite stateVector based on an Ephemeris
{

#ifdef DEBUG
cout << "In Function:  getEphemeris" << endl;
#endif

	const int DIM = 3;
	const int N_INTERP = 5;
        const double delta = 0.000001; // 1 microsec
	int i, j;
	int arrindex; 
	double timeOffset;
	double time2Interp;
	double *posArr;
	double *velArr;

	ephStart = etime[0];
	numPts = etime.size();
	ephEnd = etime[numPts-1];
	timeStep = etime[1]-etime[0];
	coordSys = ephem.coordSys;	

	//Check TimeStep and Ephemeris Size match
	double ephTS = (ephem.ephEnd - ephem.ephStart) * 86400.0/ (ephem.numPts - 1);
	if(fabs(ephTS -ephem.timeStep) > delta) {
		cerr << "ephTS( " << ephTS << " ) does not equal ephem.timeStep( " << ephem.timeStep << " )." << endl;
		cerr << "Check Ephemeris File." << endl;
		return;
	}


	i = ephem.numPts;
	//Load and Fill ephemeris array for Interpolation
	posArr = new double[i*3];
	velArr = new double[i*3];
	for(i=0; i<ephem.numPts; i++)
		for(j=0; j<3; j++)
		{
			posArr[i*3+j] = ephem.stateVector[i].position[j];
			velArr[i*3+j] = ephem.stateVector[i].velocity[j];
		}

	//Loop through each etime and interpolate ephem onto stateVector
	for(i=0; i<etime.size(); i++){
		stateVector[i].time = etime[i]; 
		calculateTime(i); 
		time2Interp = (etime[i] - ephem.ephStart)*86400.0;
		arrindex =  time2Interp / ephTS + N_INTERP/2;
		timeOffset = time2Interp - double(arrindex) * ephTS;
                interpolate(posArr, velArr ,arrindex, DIM, ephem.numPts, ephTS, N_INTERP, timeOffset, stateVector[i].position, stateVector[i].velocity);
/* // Interpolation Details
		cout << i << "  Time:  " << stateVector[i].time << "  T2I:  " << time2Interp;
		cout << "  arrind:  " << arrindex << " tOffset: " << timeOffset << endl;
		cout << i << "  Pos:  " << stateVector[i].position[0] << "   ";
		cout << stateVector[i].position[1] <<  "  " << stateVector[i].position[2] << endl;
		cout << i << "  Vel:  " << stateVector[i].velocity[0] << "   ";
		cout << stateVector[i].velocity[1] <<  "  " << stateVector[i].velocity[2] << endl;
*/ 

	}
	delete posArr;
	delete velArr;
 	return;

}


void Satellite::getTPLocation(Ephemeris& T)
//!< Gets the XYZ Tangent Point Location for each time from the its ephemeris to another ephemeris
//!< Requires an Ephemeris
{
#ifdef DEBUG
cout << "In Function:  getTPLocation" << endl;
#endif
	int i, j;
        const int DIM = 3;
        const int N_INTERP = 5;
        const double delta = 0.000001; // 1 microsec
        int arrindex;
        double timeOffset;
        double time2Interp;
        double *posArr;
        double *velArr;
	double pos_Sun[3];
	double vel_Sun[3];
	GeoCal geo;
	DVector *SC_Los_Eci, *SC_Los_Eci_UV;
	DVector *SC_Los_Ecf, *SC_Los_Ecf_UV;
	DVector *SC_Pos_Eci, *SC_Pos_Eci_Offset;
	//DVector SC_Pos_Eci_Offset(3);
	//DVector SC_Pos_Eci(3), Sun_Pos_Eci(3);
	DVector *SC_Pos_Ecf, *SC_Pos_Ecf_Offset;
	DVector *Sun_Pos_Eci, *SC_TP_Ecf;
	DVector *SC_TP_Eci;
	//DVector *SC_TP_Ecf;
	//DVector *SC_TP_Ecf;
	SC_Los_Eci = new DVector(3);
	SC_Los_Eci_UV = new DVector(3);
        SC_Pos_Eci_Offset = new DVector(3);
        SC_Pos_Ecf = new DVector(3);
        SC_Pos_Ecf_Offset = new DVector(3);
        SC_Los_Ecf = new DVector(3);
        SC_Los_Ecf_UV = new DVector(3);
	SC_TP_Ecf = new DVector(3);
	SC_TP_Eci = new DVector(3);
	Sun_Pos_Eci = new DVector(3);
	SC_Pos_Eci = new DVector(3);

        i = T.numPts;
        posArr = new double[i*3];
        velArr = new double[i*3];

        for(i=0; i<T.numPts; i++)
                for(j=0; j<3; j++)
                {
                        posArr[i*3+j] = T.stateVector[i].position[j];
                        velArr[i*3+j] = T.stateVector[i].velocity[j];
                }

	//Loop through Satellite Object Times
        for(i=0; i<numPts; i++){
		//DVector *Sun_Pos_Eci, *SC_Pos_Eci;
                //stateVector[i].time = etime[i];
                time2Interp = (stateVector[i].time - T.ephStart) * 86400.0;
                arrindex =  time2Interp / T.timeStep + N_INTERP/2;
                timeOffset = time2Interp - double(arrindex) * T.timeStep;
                //timeOffset = double(arrindex) * T.timeStep - time2Interp;
//                jdntimetest = starttime + (timestep * arrindex)/86400.0;
                //interpolate(ephem.stateVector, arrindex, DIM, asize, timestep, N_INTERP, timeoffset, stateVector.position, stateVector.velocity);
                interpolate(posArr, velArr, arrindex, DIM, T.numPts, T.timeStep, N_INTERP, timeOffset, pos_Sun, vel_Sun);

		int syear = stateVector[i].year;
		int sdoy = stateVector[i].doy;
		double ssec = stateVector[i].sec;

		Sun_Pos_Eci->init(pos_Sun, 3);
		SC_Pos_Eci->init(stateVector[i].position, 3);
		//*SC_Los_Eci = *SC_Pos_Eci - *Sun_Pos_Eci;
		*SC_Los_Eci = *Sun_Pos_Eci - *SC_Pos_Eci;
		objectDistance[i] = SC_Los_Eci->getMag();
		earthSunDistance[i] = Sun_Pos_Eci->getMag();
                //cout << "OD:  " << objectDistance[i] << "  ESD:  " << earthSunDistance[i] << endl;

		//*losVectors[i] = *SC_Los_Eci;
		*SC_Los_Eci_UV = *SC_Los_Eci;
		SC_Los_Eci_UV->unitVector();
		
		

                double sdec, srasn, gstmean;
		int notan;
		double tplatECI, tpLonECI, tpAltECI;
                //geo.cal_mean_greenwich_sidereal_time(orbitYear,orbitDay,timeSecDouble,sdec,srasn,gstmean);
                geo.cal_mean_greenwich_sidereal_time(syear,sdoy,ssec,sdec,srasn,gstmean);
		//*SC_TP_Eci = geo.calTpECI(*SC_Los_Eci_UV, *SC_Pos_Eci, gstmean, tplatECI, tpLonECI, tpAltECI, notan);
		//*SC_TP_Eci = geo.calTpECI(*SC_Los_Eci_UV, *SC_Pos_Eci, gstmean, tpLat[i], tpLon[i], tpAlt[i], notan);

		*SC_TP_Eci = geo.calTpECI(*SC_Los_Eci_UV, *SC_Pos_Eci, notan);
		geo.calCart2LLA_F1999(*SC_TP_Eci, gstmean, tpLat[i], tpLon[i], tpAlt[i], notan);
		geo.calCart2LLA_F1999(*SC_Pos_Eci, gstmean, scLat[i], scLon[i], scAlt[i], notan);

		//Debugging the Disconinuity at the Geoid`
		//cout << "(Y, D, S):  " << syear << ", " << sdoy << ", " << ssec ;
		//cout << "  (d, r, gmean):  " << sdec << ", " << srasn << ", " << gstmean ;
		//cout << "  notan:  " << notan ;
		//cout << "Lat:  " << tplatECI << "  Lon:  " << tpLonECI << "  Alt:  " << tpAltECI << endl;


		*losVectors[i] = *SC_TP_Eci - *SC_Pos_Eci;
                *losSunVectors[i] = *SC_Los_Eci;
		//losVector[i].los = SC_Los_Eci_UV;
		//*losVectors[i] = *SC_Los_Eci_UV;
		//losVectors[i]->cv[0] = SC_Los_Eci_UV->cv[0];
		//losVectors[i]->cv[1] = SC_Los_Eci_UV->cv[1];
		//losVectors[i]->cv[2] = SC_Los_Eci_UV->cv[2];
		//losVectors[i]->init(SC_Los_Eci_UV->cv, 3);
		*SC_Pos_Eci_Offset = *SC_Pos_Eci + *SC_Los_Eci_UV;
		*SC_Pos_Ecf = geo.convert_ECItoECF(syear,sdoy,ssec,*SC_Pos_Eci);
		*SC_Pos_Ecf_Offset = geo.convert_ECItoECF(syear,sdoy,ssec,*SC_Pos_Eci_Offset);
		//*SC_Los_Ecf = *SC_Pos_Ecf - *SC_Pos_Ecf_Offset;
		*SC_Los_Ecf = *SC_Pos_Ecf_Offset - *SC_Pos_Ecf;
		*SC_Los_Ecf_UV = *SC_Los_Ecf; // Currently Undefined
		SC_Los_Ecf_UV->unitVector(); // Currently Undefined
		*SC_TP_Ecf = geo.calTpECF(*SC_Los_Ecf_UV, *SC_Pos_Ecf);
		//cout << "SCTPECF:  "  << SC_TP_Ecf->cv[0] << "  I:   "  << i << endl;
		//*SC_Los_Ecf = *SC_Pos_Ecf - *SC_TP_Ecf;
		*SC_Los_Ecf = *SC_TP_Ecf - *SC_Pos_Ecf;
                scRad[i] = SC_Pos_Eci->getMag();
		//*losVectors[i] = *SC_Los_Ecf;
		Sun_Pos_Eci->clear();
		SC_Pos_Eci->clear();
		//geo.calECF2Geodetic(*SC_TP_Ecf, tpLat[i], tpLon[i], tpAlt[i]);
               //         jdntimetest = starttime + (timestep * arrindex + timeoffset)/86400;
               //         jdn2cal(jdntimetest, caltimetest);
        }
	//Interpolate Ephemeris T position
	//Calculate LOS ECI
	//
        delete posArr;
        delete velArr;
	return;
}
                                                                                                                                                             
void Satellite::getViewingAngle()
//!< Gets the Viewing Angle for each time stamp
//!< Requires LOSVector to be populated.
{
#ifdef DEBUG
cout << "In Function:  getViewingAngle" << endl;
#endif

        const double PICONST = 3.141592653589793;
	double theta;
        DVector position(3);
        DVector velocity(3);
        DVector x(3);
        DVector y(3);
        DVector z(3);
        DVector normal(3);
        //DVector *LOSUV;
	//LOSUV = new DVector(3);
	DVector LOSUV(3);
	double lossc[3];
	double hyp;

	double one, two, three;
	for(int i=0; i<numPts; i++){
		position.init(stateVector[i].position, 3);
		velocity.init(stateVector[i].velocity, 3);
		position.unitVector();
		velocity.unitVector();
		y = position ^ velocity;
		y.unitVector();
		x = y ^ position;
		x.unitVector();
		y = position ^ x;
		y.unitVector();
		LOSUV = (*losVectors[i]);
		//LOSUV = losVector[i].los;
		//LOSUV = losVectors[i];
		//LOSUV->unitVector();
		LOSUV.unitVector();
		theta = position * LOSUV;
		theta = acos(theta);
		theta = PICONST - theta;
		//cout << "Theta:  " << theta << endl;
		lossc[0] = x*LOSUV;
		lossc[1] = y*LOSUV;
		lossc[2] = position*LOSUV;
		hyp = sqrt(lossc[0]*lossc[0]+lossc[1]*lossc[1]);
		//elevation[i] = atan(lossc[2] / hyp)/RADperDEGREE; 
		elevation[i] = atan(lossc[2] / hyp); 
		position.clear();
		LOSUV.clear();
		//LOSUV->clear();
		velocity.clear();
	}
}
                                                                                                                                                           
void Satellite::calcLLA()
//!< Calculates the Geodetic Lat, Lon, and Altitude for Each
//!< Requires No Parameters at this Time
//!< IS NOT IMPLEMENTED
{

}
                                                                                                                                                             
void Satellite::calcCurvatureRadius()
//!< Calculates the Curvature Radius by Propagating Through the TP Lat, Lon, Alt Arrays
//!< Requires No Parameters at this Time
//!< IS NOT IMPLEMENTED
{
#ifdef DEBUG
cout << "In Function:  calcCurvatureRadius" << endl;
#endif

   int indCR, indCloud;
   double Rmeri, Rprime, heading, radCurvTest;
   double surfaceHeight = 0.0;
   double cloudHeight = 83.0;

   getIndex(surfaceHeight, indCR);

   calcCurvatureValues(indCR, Rprime, Rmeri);

   calcBearing(indCR, heading);

   radCurvature = Rmeri*Rprime /
                  (Rprime*cos(heading)*cos(heading)+Rmeri*sin(heading)*sin(heading));


   solarExtent = atan(SUNRAD / objectDistance[indCR]) * 2.0;
   solarRadius = earthSunDistance[indCR];
   cout.setf(ios::scientific,ios::floatfield);
   cout.precision(10);
   cout << "Bearing for Altitude (" << surfaceHeight << " Km):  " << heading << " rad  OR " << heading/RADperDEGREE << " deg" << endl;
   cout << "Rcurve:  " << radCurvature << endl;
   cout << "sol2ec:  " << solarRadius << endl;
   cout << "sol2sc  " << objectDistance[indCR] << endl;


   //Add 83Km to output
   getIndex(cloudHeight, indCloud);

   calcCurvatureValues(indCloud, Rprime, Rmeri);

   calcBearing(indCloud, heading);

   losBearing = heading;

   radCurvTest = Rmeri*Rprime /
                  (Rprime*cos(heading)*cos(heading)+Rmeri*sin(heading)*sin(heading));

   cout << "Bearing for Altitude (" << cloudHeight << " Km):  " << heading << " rad  OR " << heading/RADperDEGREE << " deg" << endl;
   cout << "Rcurve:  " << radCurvTest << endl;

   //Calculate SCRadius = SCAlt at SCLat + radCurvature
   //Initial Cut for Quick Release
   double oldRad;
   for(int i=0; i<numPts; i++){
      //scRad[i] = scAlt[i] + radCurvature;
      oldRad = scAlt[i] + radCurvature;
      scRad[i] = sqrt((losVectors[i]->getMag()*losVectors[i]->getMag())+(tpAlt[i] + radCurvature)*(tpAlt[i] + radCurvature));
      //cout << "SCRAD:  " << scRad[i] << "  " << oldRad << "  " << scRad[i] - oldRad << endl;
   }

   return;
}

void Satellite::calcBearing(int ind, double& bearing)
//!< Calculates the Bearing from SC to TP
{
#ifdef DEBUG
cout << "In Function:  calcBearing" << endl;
#endif


   double londif = scLon[ind] - tpLon[ind];
   if(londif > 180.0) londif -= 360.0;     
   while(londif < -180.0) londif += 360.0 ;
   cout << "scLat:  " << scLat[ind] << "  scLon:  " << scLon[ind] << endl;
   cout << "LonDif:  " << londif;
   double opp = -sin(londif*RADperDEGREE)*cos(tpLat[ind]*RADperDEGREE);
   cout << "  Opp:  " << opp;
   double adj = cos(scLat[ind]*RADperDEGREE)*sin(tpLat[ind]*RADperDEGREE)
                -sin(scLat[ind]*RADperDEGREE)*cos(tpLat[ind]*RADperDEGREE)*cos(londif*RADperDEGREE);
   cout << "  Adj:  " << adj << endl;
   bearing = fmod(atan2(opp,adj),2*PI);

   return;
}

void Satellite::calcCurvatureValues(int ind, double &prime, double &meridional)
//!< Calculates the Curvature Values needed
{
#ifdef DEBUG
cout << "In Function:  calcCurvatureValues" << endl;
#endif

   double ecc, amajor;
   ecc = 0.08181919084;
   amajor = 6378.137;
   double e2 = ecc*ecc;

   double sinLat = sin(tpLat[ind]*RADperDEGREE);
   double sinLat2 = pow(sinLat, 2);
   meridional = amajor*(1 - e2) / pow(1-e2*sinLat2, 1.5);
   prime = amajor /  sqrt(1-e2*sinLat2);

   cout << "TPLat:  " << tpLat[ind] << "  Rmeri: " << meridional << "  Rprime:  " << prime << endl;

  return;

}

void Satellite::getIndex(double alt, int &ind)
//!< Returns Index at input alt
{
#ifdef DEBUG
cout << "In Function:  getIndex" << endl;
#endif

   EventVar<double> CCR_TPAlt("CCR_TPAlt", "TPAlt for Curvature Use", tpAlt, numPts);

   ind = search(CCR_TPAlt, alt);
   if(ind == -1){  
      cout << "WARNING:  Unable to locate Tangent Point at Zero Km." << endl;
      cout << "          Valid Range:  [" << CCR_TPAlt[0] << ", " << CCR_TPAlt[CCR_TPAlt.size()-1] << "]" << endl;
      ind = CCR_TPAlt.size() / 2;
      cout << "          Using Tangent Point: " << CCR_TPAlt[ind] << endl;
   }

   return;

}

void Satellite::calcSCRadius()
{
   
#ifdef DEBUG
cout << "In Function:  calcSCRadius" << endl;
#endif

        const double PICONST = 3.141592653589793;
        double theta;
        DVector position(3);
        DVector velocity(3);
        DVector x(3);
        DVector y(3);
        DVector z(3);
        DVector normal(3);
        //DVector *LOSUV;
        //LOSUV = new DVector(3);
        DVector LOSUV(3);
        double lossc[3];
        double hyp;

        double one, two, three;
        for(int i=0; i<numPts; i++){
                position.init(stateVector[i].position, 3);
                velocity.init(stateVector[i].velocity, 3);
                position.unitVector();
                velocity.unitVector();
                y = position ^ velocity;
                y.unitVector();
                x = y ^ position;
                x.unitVector();
                y = position ^ x;
                y.unitVector();
                LOSUV = (*losVectors[i]);
                //LOSUV = losVector[i].los;
                //LOSUV = losVectors[i];
                //LOSUV->unitVector();
                LOSUV.unitVector();
                theta = position * LOSUV;
                theta = acos(theta);
                theta = PICONST - theta;
                //cout << "Theta:  " << theta << endl;
                lossc[0] = x*LOSUV;
                lossc[1] = y*LOSUV;
                lossc[2] = position*LOSUV;
                hyp = sqrt(lossc[0]*lossc[0]+lossc[1]*lossc[1]);
                //elevation[i] = atan(lossc[2] / hyp)/RADperDEGREE;
                elevation[i] = atan(lossc[2] / hyp);
                position.clear();
                LOSUV.clear();
                //LOSUV->clear();
                velocity.clear();
        }

return;

}



void Satellite::calcDopplerShift()
//!< Calculates Doppler Shift and stores in doppler variable
//!< Requires stateVector and losVectors populated
{
#ifdef DEBUG
cout << "In Function:  calcDopplerShift" << endl;
#endif
	for(int i=1; i<numPts; i++){
		doppler[i] = (losVectors[i]->getMag() - losVectors[i-1]->getMag()) / 
					((stateVector[i].time - stateVector[i-1].time)*86400.0);
	}
        doppler[0] = doppler[1];	

}

void Satellite::calcSinkRate()
//!< Calculates Doppler Shift and stores in doppler variable
//!< Requires stateVector and elevation variables populated
{
#ifdef DEBUG
cout << "In Function:  calcSinkRate" << endl;
#endif

   //Added to make elevation consistent with RadCurvature Method - CWB V01.02
   double ratio;
   for(int i=0; i<numPts; i++){
      ratio = (radCurvature + tpAlt[i]) / scRad[i];
      elevation[i] = acos(ratio);
   }
   //End Addition

   for(int i=1; i<numPts; i++){
      sinkRate[i] = (elevation[i] - elevation[i-1]) /
                    ((stateVector[i].time - stateVector[i-1].time)*86400.0);
   }
   sinkRate[0] = sinkRate[1];

   return;
}
                                                                                                                                                            
void Satellite::writeOutput(Event& L1, Event& Tmp)
//!< Writes the Output of All Member Variables to Event Object Tmp or Level1
//!< Requires No Parameters ... Uses Config for Output Parameters and Locations
//!< IS NOT IMPLEMENTED
{
#ifdef DEBUG
cout << "In Function:  writeOutput" << endl;
#endif


/*   EventVar <double> L1_SVPosX(numPts);
   EventVar <double> L1_SVPosY(numPts);
   EventVar <double> L1_SVPosZ(numPts);
   EventVar <double> L1_SVVelX(numPts);
   EventVar <double> L1_SVVelY(numPts);
   EventVar <double> L1_SVVelZ(numPts);
*/


   //!< Elevation Angle
   EventVar <double> L1_ViewingAngle("L1_ViewingAngle", "Elevation Angle", "Radians", elevation, numPts) ;
   //!< Tangent Point Latitude
   EventVar <double> L1_TPLat("L1_TPLat", "Tangent Point Latitude", "Degrees", tpLat, numPts);
   //!< Tangent Point Longitude
   EventVar <double> L1_TPLon("L1_TPLon", "Tangent Point Longitude", "Degrees", tpLon, numPts);
   //!< Tangent Point Altitude
   EventVar <double> L1_TPAlt("L1_TPAlt", "Tangent Point Altitude", "Km", tpAlt, numPts);
   //!< Tangent Point Curvature of Radius
   EventVar <double> L1_CurvatureRadius("L1_CurvatureRadius", "Radius of Curvature", "Km", radCurvature, 1);
   //!< Tangent Point Curvature of Radius
   EventVar <double> L1_Bearing("L1_BearingLOS", "LOS Vector Heading", "CW From N", losBearing, 1);
   //!< Doppler Effect
   EventVar <double> L1_DopplerShift("L1_DopplerShift", "Doppler Shift", "Km/s", doppler, numPts);
   //!< Elevation Angle Sink Rate 
   EventVar <double> L1_SinkRate("L1_SinkRate", "SinkRate", "Rad/s", sinkRate, numPts);
   //!< Solar Distance 
   EventVar <double> L1_SolarDistance("L1_SolarDistance", "Solar Distance From Satellite Position", "Km", objectDistance, numPts);
   //!< Solar Distance 
   EventVar <double> L1_EarthSunDistance("L1_EarthSunDistance", "Solar Distance From Earth Center", "Km", solarRadius, 1);
   //!< Spacecraft Geodetic Latitude
   EventVar <double> L1_SCLat("L1_SCLat", "Spacecraft Geodetic Latitude", "Degrees", scLat, numPts);
   //!< Spacecraft Geodetic Longitude
   EventVar <double> L1_SCLon("L1_SCLon", "Spacecraft Geodetic Longitude", "Degrees", scLon, numPts);
   //!< Spacecraft Geodetic Altitude
   EventVar <double> L1_SCAlt("L1_SCAlt", "Spacecraft Geodetic Altitude", "Km", scAlt, numPts);
   //!< Spacecraft Radius from Center of Curvature
   EventVar <double> L1_SCRad("L1_SCRad", "Spacecraft Radius from Center of Curvature", "Km", scRad, numPts);
   //!< Solar Extent
   EventVar <double> L1_SolarExtent("L1_SolarExtent", "Solar Extent", "Km", solarExtent, 1);


//EventVar<double> J("Name2",data,n);


   #ifdef MAKEPLOTS
   EventVar <double> timeArray;
   //L0.getEventVar("L0_EventTimes", timeArray);
   Tmp.getEventVar("TReg_timeJDN", timeArray);
   //Tmp.getEventVar("Detector Time", timeArray);
   timeArray -= ephStart;
   timeArray *= 86400.0;
   timeArray.setUnits("Secs");

   timeArray.setAxisRange();
   L1_TPLat.setAxisRange();
   L1_TPLon.setAxisRange();
   L1_TPAlt.setAxisRange(-100.0, 300.0);
   L1_SCLat.setAxisRange();
   L1_SCLon.setAxisRange();
   L1_SCAlt.setAxisRange(500.0, 700.0);
   L1_SCRad.setAxisRange(6900.0, 7100.0);
   L1_DopplerShift.setAxisRange(-10.0, 10.0);
   L1_ViewingAngle.setAxisRange();
   L1_SinkRate.setAxisRange();
   L1_SolarDistance.setAxisRange();
   L1_TPLat.Plot("Time Registration","TPLat ",timeArray);
   L1_TPLon.Plot("Time Registration","TPLon ",timeArray);
   L1_TPAlt.Plot("Time Registration","TPAlt ",timeArray);
   L1_SCLat.Plot("Time Registration","SCLat ",timeArray);
   L1_SCLon.Plot("Time Registration","SCLon ",timeArray);
   L1_SCAlt.Plot("Time Registration","SCAlt ",timeArray);
   L1_SCRad.Plot("Time Registration","SCRad ",timeArray);
   L1_ViewingAngle.Plot("Time Registration","Elevation Angle ",timeArray);
   L1_TPAlt.Plot("Time Registration","Doppler Shift", L1_DopplerShift);
   L1_TPAlt.Plot("Time Registration","Sink Rate", L1_SinkRate);
   L1_TPAlt.Plot("Time Registration","Solar Distance", L1_SolarDistance);
   #endif
   #ifdef OUTPUTLLA

   cout << "Radius of Curvature:  "  << L1_CurvatureRadius[0] << endl;
   cout << "LOS Bearing  "  << L1_Bearing[0] << endl;
   cout << "Solar Extent:  "  << L1_SolarExtent[0] << endl;
   cout << "Time\tElevation\tLat\tLon\tAlt\tSCRad" << endl;
   for(int i=0; i<numPts; i++){
	cout << timeArray[i] << "\t";
	cout << L1_ViewingAngle[i] << "\t";
	cout << L1_TPLat[i] << "\t";
	cout << L1_TPLon[i] << "\t";
	cout << L1_TPAlt[i] << "\t";
	cout << L1_SCRad[i] << endl;
   }
   #endif

   L1.addEventVar(L1_TPLat);
   L1.addEventVar(L1_TPLon);
   L1.addEventVar(L1_TPAlt);
   L1.addEventVar(L1_DopplerShift);
   L1.addEventVar(L1_ViewingAngle);
   L1.addEventVar(L1_SinkRate);
   L1.addEventVar(L1_SolarDistance);
   L1.addEventVar(L1_EarthSunDistance);
   L1.addEventVar(L1_CurvatureRadius);
   L1.addEventVar(L1_Bearing);
   L1.addEventVar(L1_SCLat);
   L1.addEventVar(L1_SCLon);
   L1.addEventVar(L1_SCAlt);
   L1.addEventVar(L1_SCRad);
   L1.addEventVar(L1_SolarExtent);
   
   
   return;

}