//
//  $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:	ZPTProfile.cpp
//
// Author:	John Burton	<john@gats-inc.com>
//
// Date:  	Thu May  4 17:01:31 2006
//
//-----------------------------------------------------------------------
//
// Modification History:
//
//     $Log$
//
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

#include "ConfigFile.h"
#include "EventVar.h"
#include "ZPTProfile.h"
#include "Assemble_ZPT.h"
#include "SolarFlux.h"
#include "EphemData.h"
#include "EventPopul.hpp"
#include "UTC_Conv.h"
#include "GATS_Exception.h"


#include "GATS_Utilities.hpp"


#include <iostream>
#include <cassert>
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//



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



//
//-----------------------------------------------------------------------
// Utility Routines:
//-----------------------------------------------------------------------
//
void ConvValArrayToDouble(valarray<float> const &input, valarray<double> &output);
void Check_Fix_NCEP( valarray<double>& z_ncep, valarray<double>& p_ncep, valarray<double>& t_ncep);

int ZPTProfile(Event& L0, Event& L1, Event& Tmp, Event& SD, ConfigFile& cf)
{

const string MODULE_NAME = "ZPTProfile";

EventVar <float> workFltArray; 
//EventVar <int>    workIntArray; 
EventVar <double> Z_Var;  
EventVar <float> flt_Var;  
EventVar <double> P_Var;  
EventVar <double> T_Var;  
EventVar <double> Lat_Var;  
EventVar <double> Lon_Var;  
//EventVar <long> KP_Array;  
EventVar <int> KP_Array;  

//EventVar <double> workFltArray;  

SolarFlux flux; 
EphemData ephemstd;
Assemble_ZPT assembledProf;

valarray<double> z_merged; 
valarray<double> p_merged; 
valarray<double> t_merged;
valarray<double> work;

int retCode = 0;

try {

    Tmp.getEventVar("NCEP Ztp", Z_Var);
    Tmp.getEventVar("NCEP Pressure", flt_Var); 
    ConvValArrayToDouble(flt_Var, P_Var);
    Tmp.getEventVar("NCEP Temperature", flt_Var); 
    ConvValArrayToDouble(flt_Var, T_Var);

    valarray<double> z_ncep = Z_Var; 
    valarray<double> p_ncep = P_Var; 
    valarray<double> t_ncep = T_Var;

Check_Fix_NCEP(z_ncep, p_ncep, t_ncep) ;

    L1.getEventVar("L1_TPAlt", Z_Var);
    L1.getEventVar("L1_TPLat", Lat_Var);
    L1.getEventVar("L1_TPLon", Lon_Var);

    valarray<double> ephem_alt = Z_Var;
    valarray<double> ephem_lat = Lat_Var;
    valarray<double> ephem_lon = Lon_Var;

    valarray<int> kp(1, KP_ARRAY_SIZE); // hardwire this to 1 values until further along in integration 
    
    Tmp.getEventVar("SGI Kp Index", workFltArray);
    assert(workFltArray.size() == kp.size());
    for(int i = 0; i < workFltArray.size(); i++)
        { kp[i] = (int)workFltArray[i]; }    
    
    Tmp.getEventVar("SGI Ap Index", workFltArray); 
    float flt_temp = workFltArray[0];    
    long apvg = flt_temp; 

    Tmp.getEventVar("SGI Solar Flux", workFltArray); 
    float f10_7 = workFltArray[0];

    float f_90 = 150.0; // not part of the Aux data. Hardwired to SABER processing's value
    
    time_t timeIn = static_cast<long>(L1.getEventStartTime());
    UTC_Conv  startTime(timeIn);
    
    int date = (startTime.Year()%100)*1000 + startTime.Doy(); 
    double sec = (startTime.Hour()*60. + startTime.Min())*60. + startTime.Sec();

    flux = SolarFlux(kp, apvg, f10_7, f_90);
    ephemstd = EphemData(ephem_alt, ephem_lat, ephem_lon);
//std::cout << " befor assemble " << z_ncep.size() << std::endl;
    assembledProf = Assemble_ZPT (date, sec, flux, ephemstd, z_ncep, 
                                                        p_ncep, t_ncep);
//std::cout << "after assemble " << std::endl;
    assembledProf.getAltitude(z_merged);

    EventPopul(z_merged, "Z_Merged", L1);

    assembledProf.getPressure(p_merged);
    EventPopul(p_merged, "P_Merged", L1);

    assembledProf.getTemperature(t_merged);
    EventPopul(t_merged, "T_Merged", L1);
/*
for (int i = 0; i< z_merged.size() ; ++i ) {


std::cout << p_merged[i] << "  " << z_merged[i] << std::endl;

}
exit(23) ;
*/

    assembledProf.getXRatioAlt(work);
    EventPopul(work, "MSIS_Alt", L1);
    
    assembledProf.getxRatio_O(work);
    EventPopul(work, "MSIS_O", L1);
    
    assembledProf.getxRatio_O2(work);
    EventPopul(work, "MSIS_O2", L1);
    
    assembledProf.getxRatio_N2(work);
    EventPopul(work, "MSIS_N2", L1);
    
    } catch (GATS_Exception gex) {
    
        cerr << "GATS Exception Caught in " << MODULE_NAME << ": " << gex.what() << endl;
        retCode = -300;
        
    } catch (exception ex) {

        cerr << "Low Level Exception Caught in " << MODULE_NAME << ": " << ex.what() << endl;
        retCode = -999;    
    }

//  std::cerr << "In ZPTProfile \n";
  
  return retCode;
}

void ConvValArrayToDouble(valarray<float> const &input, valarray<double> &output)
{
output.resize(input.size());

for(unsigned long i = 0; i < input.size(); i++)
    output[i] = input[i];    

}


void Check_Fix_NCEP( valarray<double>& z_ncep, valarray<double>& p_ncep, valarray<double>& t_ncep)
{
using namespace GATS_Utilities;

   assert(z_ncep.size() == p_ncep.size() &&  z_ncep.size() == t_ncep.size() );

   int N = (int)z_ncep.size();


// find negative hts at the top of the profile  (this can happen for missing data) and set a new N
//std::cout << " size " << N << " " << z_ncep[0] << "  " << z_ncep[3] << std::endl ;

//std::cout << monotonic_if( &z_ncep[0], &z_ncep[10],std::less<double>() ) << std::endl;

   while( N>2 && !monotonic_if( &z_ncep[0], &z_ncep[N],std::less<double>() )  )  N--;

   assert(N > 2);

// check that pressures are all greater than 1.e-3
   assert(&p_ncep[N] ==  find_if(&p_ncep[0], &p_ncep[N], std::bind2nd(std::less<double>(), 1.e-3) ) );

// check that pressure are monotonically decreasing
   assert (monotonic_if( &p_ncep[0], &p_ncep[N],std::greater<double>() ) );

// check that temperatures are all > 100K
   assert(&t_ncep[N] ==  find_if(&t_ncep[0], &t_ncep[N], std::bind2nd(std::less<double>(), 100.) ) );

// now here's the fix if N < original size (ie. missing data on top) we trim that off  and
// reset the val arrays 

  	if(z_ncep.size() != (int)N ) {
        valarray<double> sht(&z_ncep[0], N); 
		z_ncep.resize(N);
        z_ncep = sht;
        sht = valarray<double>(&t_ncep[0], N);
        t_ncep.resize(N);
        t_ncep = sht;
        sht = valarray<double>(&p_ncep[0], N);
        p_ncep.resize(N);
		p_ncep = sht;
   }
//  std::cout << " check ncep " << z_ncep.size() << "  " << p_ncep.size() << "  " << t_ncep.size() <<  std::endl;


}