//
//  $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.
//
//-----------------------------------------------------------------------
//
// Module:	AuxDataInput.cpp
//
// Author:	John Burton	<john@gats-inc.com>
//
// Date:  	Thu May  4 17:01:31 2006
//
//-----------------------------------------------------------------------
//
// Modification History:
//
//    
// Description:
//	Retrieves and Stores SGI and NCEP data in EventVars
//      NetCDF Files from TIMED have folliwing info:  (source: http://www.timed.jhuapl.edu/mdc/data/)
//	   SGI Data:
//	      Daily Values for Ap, 8 Kps  (source: http://sec.noaa.gov/ftpdir/latest/DGD.txt)
//                     and observed F10.7 (source: http://www.drao.nrc.ca/icarus/www/current.txt)
//         NCEP Data:
//            PressureLevel, GeoPotentialHeight, Temperature at 
//            Lat on Integer Range: [-90, 90]
//            Lon on Integer Range: [0, 359]
//            Pressure on Range: { 1000, 850, 700, 500, 400, 300, 250, 200, 
//                                 150, 100, 70, 50, 30, 10, 5, 2, 1, 0.4 }
//
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

//#include "ConfigFile.h"
//#include "EventVar.h"
#include "AuxDataInput.h"
//#include "NMCdata.h"
#include "GeneralFunctions.h"
#include "SGI_FileReader.h"
#include "NMC_FileReader.h"
#include "NMCdata.h"
#include "GATS_Utilities.hpp"
#include <iostream>
#include <string>
#include <exception>
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//
//#define TPEPSILON 0.1
//#define DOUBLE_EQ(x,v) (((v-TPEPSILON)<x)&&((v+TPEPSILON)>x))

using GATS_Utilities::ConvertToString;
//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//



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

int AuxDataInput(Event& L0, Event& L1, Event& Tmp, Event& SD, ConfigFile& cf)
{
   int i;
   int nmcStatus;
   //double nmcTP;
   char* auxDir;
   int dateCal[7];
   int sgiDateCal[7];
   int nmcDateCal[7];
   //char sgiFile1[35];
   //char sgiFile2[35];
   char *sgiFile1;
   char *sgiFile2;
   char nmcVersion[]="01";
   char nmcHour[]="12";
   char nmcFileName[140];
   char sgiFileName1[140];
   char sgiFileName2[140];
   int  sgiFlag = 0;
   int  nmcFlag = 0;
   int  oLen;
   int  orbitDate;
   int  sgiOrbitDate;
   int  nmcOrbitDate;
   int  sgiAge, nmcAge;
   double orbitDateJDN;
   double sgiOrbitDateJDN;
   double nmcOrbitDateJDN;
   char orbitDateTxt[8];
   char nmcOrbitDateTxt[7];
   char endDate_SGI1[8];
   char begDate_SGI2[8];
   char endDate_SGI2[8];
   double eventStartTime;
 
   EventVar <double> ADI_timeArray;
   EventVar <double> ADI_Lat;
   EventVar <double> ADI_Lon;
   EventVar <double> ADI_Alt;
   EventVar <int> SGI_Age("SGI_Age", "Age of Solar Geomagnetic Indices File Used", "Days", -999, 1) ;
   EventVar <int> NCEP_Age("NCEP_Age", "Age of NCEP Data File Used", "Days", -999, 1) ;

   SGI_FileReader *sgi1;
   SGI_FileReader *sgi2;
   NMC_FileReader *nmc;
   
   string strErr;
   
   //Get Config File Parameters
   //NetCDF Naming Conventions:
   //   SGI File: MDC__SYYYDOY_EYYYDOY_VV.ind
   //             where SYYYDOY is Startdate
   //             EYYYDOY is the Enddate and
   //             VV is the Version
   //   Example:  MDC__2005121_2005169_01.ind
   //   
   //   NCEP File:MDC__YYYYDOY12_VV.cpc
   //             where YYYYDOY is date and 
   //             VV is the Version 
   //   Example:  MDC__200219412_01.cpc


   const int PLOTOUTPUT = cf.GetInt("AuxDataInput", "Plot_Output");
   const int LOGOUTPUT = cf.GetInt("AuxDataInput", "Log_Output");
   const double nmcTP = cf.GetReal("AuxDataInput", "NCEP_TangentAlt");
   const double maxDelta = cf.GetReal("AuxDataInput", "MaxDays2UseOldFiles");

   auxDir = cf.GetStr("AuxDataInput", "Auxilliary_Directory");
   sgiFile1 = cf.GetStr("AuxDataInput", "SGI_LongFile");
   sgiFile2 = cf.GetStr("AuxDataInput", "SGI_ShortFile");
   //nmcTP = cf.GetReal("AuxDataInput", "NCEP_TangentAlt");

   eventStartTime = L0.getEventStartTime();
   //Obtain Inputs from Event Classes
   //Get Event Variables 
   L0.getEventVar("L0_DetectorTimes", ADI_timeArray);
   L1.getEventVar("L1_TPLat", ADI_Lat);
   L1.getEventVar("L1_TPLon", ADI_Lon);
   L1.getEventVar("L1_TPAlt", ADI_Alt);

   //InitializeEventVars(dataRegEVS);

   //Get SGI File Dates
   for (i=0;i<7;i++){
       endDate_SGI1[i]=sgiFile1[13+i];
       begDate_SGI2[i]=sgiFile2[5+i];
       endDate_SGI2[i]=sgiFile2[13+i];
   }
       endDate_SGI1[7]='\0';
       begDate_SGI2[7]='\0';
       endDate_SGI2[7]='\0';

   orbitDateJDN =// ADI_timeArray[0]+
       UNIX_CONST+eventStartTime/86400.0;
   nmcOrbitDateJDN = orbitDateJDN;
   jdn2cal(orbitDateJDN, dateCal);
   orbitDate = cal2doy(dateCal);
   nmcOrbitDate = orbitDate;
   oLen = sprintf(orbitDateTxt,"%u",orbitDate);

   sprintf(nmcFileName,"%sMDC__%s%s_%s.cpc",auxDir,orbitDateTxt,nmcHour,nmcVersion);
   sprintf(sgiFileName1,"%s%s",auxDir,sgiFile1);
   sprintf(sgiFileName2,"%s%s",auxDir,sgiFile2);
   nmc = new NMC_FileReader();

  // nmc = new NMC_FileReader(nmcFileName);
   nmcStatus = 0;
   //cout << "NMCJDN:  " << nmcOrbitDateJDN << endl;
   //cout << "ORBJDN:  " << orbitDateJDN << endl;
   //cout << "ORBDAT:  " << orbitDate << endl;

   //while(!nmcStatus && orbitDateJDN < nmcOrbitDateJDN + 10.0)
   while(!nmc->getStatus() && orbitDateJDN < nmcOrbitDateJDN + maxDelta)
   {
      try {
         //nmc = new NMC_FileReader(nmcFileName);
         nmc->SetUp(nmcFileName);
      }
      catch (const char *message) {
        cout << "***File not found.  Will Try " << nmcOrbitDateJDN - orbitDateJDN + maxDelta << " more days." << endl;
	cout <<  message << endl;
	strErr = string("*** NMC File not found.  Will Try ") + 
			ConvertToString(nmcOrbitDateJDN - orbitDateJDN + maxDelta) + 
			string(" more days.");
	L1.addLogEntry(strErr);
	L1.addLogEntry(string(message));
	
        //cout << "NMCJDN:  " << nmcOrbitDateJDN << endl;
	nmcOrbitDateJDN-=1.0;
        //cout << nmcOrbitDateJDN << endl;
	jdn2cal(nmcOrbitDateJDN, nmcDateCal);
        //cout << nmcDateCal[2] << endl;
        nmcOrbitDate = cal2doy(nmcDateCal);
        //cout << nmcOrbitDate << endl;
        sprintf(nmcOrbitDateTxt,"%d",nmcOrbitDate);
        sprintf(nmcFileName,"%sMDC__%s%s_%s.cpc",auxDir,nmcOrbitDateTxt,nmcHour,nmcVersion);
        //cout << nmcFileName << endl;
      }
        nmcFlag = int (orbitDateJDN - nmcOrbitDateJDN + 1e-15);
        //cout << nmcFlag << endl;
        nmcStatus = nmc->getStatus();
        //sleep(2);
        //cout << nmcStatus << endl;
        //cout << "END LOOP" << endl;
   }

   if(!nmcStatus) throw runtime_error("*** NMC File not found within MaxDays Time Frame.\n");
   else {
      NCEP_Age[0] = nmcFlag;
      cout << "Found NCEP File." << endl;
      cout << "OrbitDate:  " << orbitDate << endl;
      cout << "SGIFile1:  " << sgiFileName1 << endl;
      cout << "SGIFile2:  " << sgiFileName2 << endl;
      cout << "NMCFile:  " << nmcFileName << endl << endl;
   }
   
   L1.addLogEntry(string("OrbitDate: ") + ConvertToString(orbitDate));
   L1.addLogEntry(string("SGIFile1:  ") + string(sgiFileName1));
   L1.addLogEntry(string("SGIFile2:  ") + string(sgiFileName2));
   L1.addLogEntry(string("NMCFile:   ") + string(nmcFileName));
   
   //SGI_Record_Buffer sgi_record_buffer1;
   //SGI_Record_Buffer sgi_record_buffer2;
   //sgi->Read_SGI_RecordsToBuffer(sgi_record_buffer);
   sgi2 = new SGI_FileReader(sgiFileName2);
   sgi2->Read_SGI_RecordsToBuffer();
   sgi1 = new SGI_FileReader(sgiFileName1);
   sgi1->Read_SGI_RecordsToBuffer();

   float Kp[32];
   float Ap;
   float sflux;

   //Check that 2nd file contains the orbit date and load values
   //Otherwise load last available days values
   if(sgi2->getEndJDN() >= orbitDateJDN && sgi2->getBegJDN() <= orbitDateJDN){
      cout << "   Using SGIFile2" << endl;
      sgiFlag = 0;
      sgiOrbitDateJDN = orbitDateJDN;
      for(int i=3; i>=0; i--){
         jdn2cal(sgiOrbitDateJDN, sgiDateCal);
         sgiOrbitDate = cal2doy(sgiDateCal);
         if(i==3){
            Ap = sgi2->record_buffer[sgiOrbitDate].Ap;
            sflux = sgi2->record_buffer[sgiOrbitDate].flux10_7;
         }
         if(sgiOrbitDateJDN > sgi2->getBegJDN())
            sgi2->getKpValues(Kp, sgiOrbitDate, i*8);
         else sgi1->getKpValues(Kp, sgiOrbitDate, i*8); 
         sgiOrbitDateJDN -= 1.0;
      }
   }
   else if(sgi1->getEndJDN() >= orbitDateJDN){
      cout << "   Using SGIFile1" << endl;
      sgiFlag=0;
      sgiOrbitDateJDN = orbitDateJDN;
      for(int i=3; i>=0; i--){
         jdn2cal(sgiOrbitDateJDN, sgiDateCal);
         sgiOrbitDate = cal2doy(sgiDateCal);
         if(i==3){
            Ap = sgi1->record_buffer[sgiOrbitDate].Ap;
            sflux = sgi1->record_buffer[sgiOrbitDate].flux10_7;
         }
         sgi1->getKpValues(Kp, sgiOrbitDate, i*8);
         sgiOrbitDateJDN -= 1.0;
      }
   }
   else {

      if((orbitDateJDN-sgi2->getEndJDN()) > maxDelta){ 
         cout << "SGI File Eclipsed MaxDelta of " <<  maxDelta << " Days." << endl;
         cout << orbitDateJDN-sgi2->getEndJDN() << " Days. " << endl;
         throw runtime_error("SGI File more than MaxDelta days old.  Edit the INI file with current SGI File.");
      }

      cout << "   Using SGIFile2:  " << orbitDateJDN-sgi2->getEndJDN() << " days past the last day." << endl;
      sgiFlag = int (orbitDateJDN - sgi2->getEndJDN() + 1.e-15);
      sgi2->getKpValues(Kp);
      Ap = sgi2->record_buffer[sgi2->date_enddoy].Ap;
      sflux = sgi2->record_buffer[sgi2->date_enddoy].flux10_7;
   }

   SGI_Age[0] = sgiFlag;

   delete sgi1;
   delete sgi2;

   EventVar <int> SGI_Date("SGI Date", orbitDate, 1);
   EventVar <int> SGI_Flag("SGI Flag", sgiFlag, 1);
   //EventVar <float> SGI_Ap("SGI Ap Index", sgi_record_buffer[orbitDate].Ap, 1);
   EventVar <float> SGI_Ap("SGI Ap Index", Ap, 1);
   //EventVar <float> SGI_Kp("SGI Kp Index", sgi_record_buffer[orbitDate].Kp[kpIndex], 32);
   EventVar <float> SGI_Kp("SGI Kp Index", Kp, 32);
   //EventVar <float> SGI_Flux107("SGI Solar Flux", sgi_record_buffer[orbitDate].flux10_7, 1);
   EventVar <float> SGI_Flux107("SGI Solar Flux", sflux, 1);

   if(LOGOUTPUT){
      cout << "SGI DATA" << endl;
      cout << "   SGIAp:  " << SGI_Ap[0] << endl;
      for(int j=0; j<32; j++)
         cout << "   SGIKp[" << j << "]:  " << SGI_Kp[j] << endl;
      cout << "   SGIFlux:  " << SGI_Flux107[0] << endl << endl;
   }
    float lat, lon;
    int nmci;
    //cout << "ADI_Alt: " << ADI_Alt[0] << "   " << ADI_Alt[ADI_Alt.size()-1] << endl;
    nmci = search(ADI_Alt, nmcTP);
    //cout << "Here." << endl;
    if(nmci == -1) {
	cout << "WARNING!!!  AuxDataIn_Alt Range: " << ADI_Alt[0] << "   " << ADI_Alt[ADI_Alt.size()-1] << endl;
        cout << "Cannot Use Nominal NMC Altitude of : "  << nmcTP << endl;
        nmci = ADI_Alt.size() / 2;
        cout << "Using Lat/Lon of 0/0" << endl;
        lat = 0.0;
	lon = 0.0;
    }
    lat = ADI_Lat[nmci];
    lon = ADI_Lon[nmci];
   float nmcPress[NMCLEVELdim],nmcGht[NMCLEVELdim],nmcTemp[NMCLEVELdim];
   double nmcZTP[NMCLEVELdim];

  if(nmc->getStatus())
    nmc->Read_netCDF_NMC_data(lat,lon,nmcPress,nmcTemp, nmcGht, nmcZTP);
  else {
    nmcFlag = -999.0;
    lat = -999.0;
    lon = -999.0;
    for(int j=0;j<NMCLEVELdim;j++){
       nmcPress[j] = -999.0;
       nmcTemp[j] = -999.0;
       nmcGht[j] = -999.0;
       nmcZTP[j] = -999.0;
    }
  }
  
  delete nmc;

    EventVar <double> NCEP_Flag("NCEP Flag", nmcFlag, 1);
    EventVar <double> NCEP_Lat("NCEP Latitude", lat, 1);
    EventVar <double> NCEP_Lon("NCEP Longitude", lon, 1);
    EventVar <float> NCEP_Press("NCEP Pressure", nmcPress, NMCLEVELdim);
    EventVar <float> NCEP_Temp("NCEP Temperature", nmcTemp, NMCLEVELdim);
    EventVar <float> NCEP_GHT("NCEP GeoPotential Height", nmcGht, NMCLEVELdim);
    EventVar <double> NCEP_ZTP("NCEP Ztp", nmcZTP, NMCLEVELdim);


    //Convert units to Kms
    NCEP_GHT /= 1000.0;
    NCEP_ZTP /= 1000.0;

   if(LOGOUTPUT){
      cout <<  "NCEP Data:  " << NCEP_Flag[0] << endl;
      cout << "   ORBDate:  " << orbitDate << endl;
      cout << "   NMCDate:  " << nmcOrbitDate << endl;
     
      cout <<  "   Latitude:  " << NCEP_Lat[0] << "\tLongitude:  " << NCEP_Lon[0] << endl;
      cout <<  "   Pressure\tTemp\tGeoHt\tZtp" << endl;
      for (int j=0;j<NMCLEVELdim;j++)
      {
             //record.pressure_nmc[j] = nmcPress[j];
             //record.temperature_nmc[j]=nmcTemp[j];
             //record.geopotential_height_nmc[j] = nmcGht[j]/1000.; //convert to km
             cout <<  NCEP_Press[j] << "\t" << NCEP_Temp[j] << "\t" << NCEP_GHT[j] << "\t" << NCEP_ZTP[j] << endl;
      }
      cout << endl;
   }

if(PLOTOUTPUT == 1){ 
    char Subtitle[50];
    lat += 0.5;
    lon += 0.5;
    lat = double( int(lat) );
    lon = double( int(lon) );
    sprintf(Subtitle, "NCEP T/P for Orbit %7i Day %7i Latitude: %5.1f  Longitude: %5.1f  ", orbitDate, nmcOrbitDate, lat, lon);
    NCEP_ZTP.setAxisRange(0, 60000);
    NCEP_ZTP.setDescription("TP Altitude");
    NCEP_ZTP.setUnits("Km");
    NCEP_Temp.setAxisRange(200, 300);
    NCEP_Temp.setDescription("Temperature");
    NCEP_Temp.setUnits("K");
    NCEP_Press.setAxisRange(0, 1000);
    NCEP_Press.setDescription("Pressure");
    NCEP_Press.setUnits("mbar");
    NCEP_Press.Plot("SOFIE AuxDataInput",Subtitle,NCEP_Temp);
    sprintf(Subtitle, "NCEP T/Z for Latitude: %5.1f  Longitude: %5.1f  ", lat, lon);
    //NCEP_ZTP.PlotVar("SOFIE AuxDataInput",Subtitle,NCEP_Temp);
}

    Tmp.addEventVar(NCEP_Lat);
    Tmp.addEventVar(NCEP_Lon);
    Tmp.addEventVar(NCEP_Press);
    Tmp.addEventVar(NCEP_Temp);
    Tmp.addEventVar(NCEP_ZTP);

    Tmp.addEventVar(SGI_Date);
    Tmp.addEventVar(SGI_Ap);
    Tmp.addEventVar(SGI_Kp);
    Tmp.addEventVar(SGI_Flux107);

    L1.addEventVar(NCEP_Age);
    L1.addEventVar(SGI_Age);

    cout << "NCEP_AGE SGI_AGE" << endl;
    cout << NCEP_Age[0] << "         " << SGI_Age[0] << endl;

    std::cerr << "In AuxDataInput \n";

    return 0;
}