#include "Level1Data.h"
#include "F77Radtran.h"
#include "GATS_Utilities.hpp"

#include <iostream>
#include <vector>
#include <stdio.h>
//#include <cstring>
//#include <cmath>

using namespace GATS_Utilities;


extern "C" void get_msis_data ( int* usemsis, int* usesaber, int* O_x2, int* naltsab, float* saberalt, float **saberNH, float **saberSH, float* zr,  float **q,  int* nr, int* igs,    int*  ngg)
{

	static int c__1 = 1;

	std::string dt = gLevel1Data.StartTime();
	struct tm timestruct={0};
	
	timestruct.tm_year = ConvertFromString<int>( dt.substr(0,4), std::dec) - 1900;
	timestruct.tm_mon = ConvertFromString<int>( dt.substr(5,2), std::dec) - 1;
	timestruct.tm_mday = ConvertFromString<int>( dt.substr(8,2), std::dec) ;

	mktime(&timestruct) ;

        float lat_temp = gLevel1Data.Lat83();

	std::cout << "Day of year: " << timestruct.tm_yday <<": Latitude: " <<lat_temp << std::endl;

	if( *usemsis == 0) return;
std::cout << "Replace O, O2 and N2 with MSIS Data " << std::endl;	
std::cout << *usemsis << " " << *usesaber << std::endl;	
	std::vector<float> R,  msisAlt = ::gLevel1Data.msisAlt();
	assert (! msisAlt.empty() );
	int n_msis = (int)(msisAlt.size());

	int* ik = std::find(igs,  igs+ *ngg , 700 );
	int j = std::distance( igs, ik);
	if(j < *ngg) {
		R = ::gLevel1Data.msisO2();
		::kp_intprof_noxtrp__( nr, &n_msis, &c__1, nr, &msisAlt[0], zr, 
			&R[0], &q[j][0]  );
	}
	ik = std::find(igs,  igs+ *ngg , 2200 );
	j = std::distance( igs, ik);
	if(j < *ngg) {
		R = ::gLevel1Data.msisN2();
		::kp_intprof_noxtrp__( nr, &n_msis, &c__1, nr, &msisAlt[0], zr, 
			&R[0], &q[j][0]  );
	}

	ik = std::find(igs,  igs+ *ngg , 3400 );
	j = std::distance( igs, ik);
std::cout << "Made it to O replacement " << std::endl;	
	if(j < *ngg)
	 {
         if( *usesaber != 0)
// Use SABER O profile in the analysis
	  {
std::cout << "Replace O with SABER Data " << std::endl;
	   std::vector<float> alt(*naltsab);
 	   std::vector<float> R(*naltsab);
	   int iday=timestruct.tm_yday;	
           if(lat_temp > 0.0) {
             for(int i = 0; i< *naltsab; ++i) {
               alt.at(i)=saberalt[i];
               R.at(i)=saberNH[iday][i];
             }
           }
           else {
             for(int i = 0; i< *naltsab; ++i) {
               alt.at(i)=saberalt[i];
               R.at(i)=saberSH[iday][i];
             }
           }
// Reverse R and alt and replace O in q array with R
           std::cout << " 1st Altitude, O from SABER" << alt[0] << " " << R[0] << std::endl;	
           std::cout << "last Altitude, O from SABER" << alt[*naltsab-1] << " " << R[*naltsab-1] << std::endl;	
           std::reverse(R.begin(), R.end());	
           std::reverse(alt.begin(), alt.end());
           std::cout << " 1st Altitude, O from SABER" << alt[0] << " " << R[0] << std::endl;	
           std::cout << "last Altitude, O from SABER" << alt[*naltsab-1] << " " << R[*naltsab-1] << std::endl;	
// If O_x2 then increase O by factor of 2
                if( *O_x2 != 0)
		{
                  std::transform(R.begin(), R.end(), R.begin(),
                  std::bind2nd( std::multiplies<float>(), 2.0  ));
                }
// OK replace values in q with R	
	   std::replace_if(R.begin() ,R.end() ,std::bind2nd(std::less<float>(),1.e-9), 1.e-9) ;
	   ::kp_intprof_noxtrp__( nr, naltsab, &c__1, nr, &alt[0], zr, 
			&R[0], &q[j][0]  );
	  }
	  else
	  {
		R = ::gLevel1Data.msisO();
// If O_x2 then increase O by factor of 2
                if( *O_x2 != 0)
		{
                  std::transform(R.begin(), R.end(), R.begin(),
                  std::bind2nd( std::multiplies<float>(), 2.0  ));
                }
		std::replace_if(R.begin() ,R.end() ,std::bind2nd(std::less<float>(),1.e-9), 1.e-9) ;
		::kp_intprof_noxtrp__( nr, &n_msis, &c__1, nr, &msisAlt[0], zr, 
			&R[0], &q[j][0]  );
	  }
	}
// comment below

/*   for (j=0;j< *ngg;++j) std::cout << igs[j] << "  " ;
    std::cout << std::endl;
	for(int i = 0; i< *nr; ++i) {

		std::cout << zr[i] ;
        for (j=0;j< *ngg;++j) std::cout << " " << q[j][i] ;
       std::cout << std::endl;
	}
	exit(0);
*/
}