#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, 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)
// Read in a SABER O profile and use that in the analysis
	  {
std::cout << "Replace O with SABER Data " << std::endl;	
	   std::vector<float> alt(100);
	   std::vector<float> R(100);
	   float dum1;
	   float dum2;
	   FILE *infile;
std::cout << "Opening SABER_O.txt " << std::endl;	
	   infile = fopen("SABER_O.txt","r");
	   if( infile == NULL)
	   {
	      fprintf(stderr, "Can't open SABER_O.txt!\n");
	      exit(1);
	   }
	   int num=0;
std::cout << num << std::endl;	
	   while (fscanf(infile, "%f %e", &dum1, &dum2) ==2)
	   {
std::cout << num << std::endl;	
std::cout << dum1 << " " << dum2 << std::endl;	
	     alt.at(num) = dum1;
	     R.at(num) = dum2;
	     num=num+1;
	   }
std::cout << num << std::endl;
	   fclose(infile);	  
std::cout << "Closed SABER_O.txt " << std::endl;
           alt.resize(num);	
           R.resize(num);	
// Reverse R and alt and replace O in q array with R
        std::cout << alt[0] << std::endl;	
        std::cout << alt[num-1] << std::endl;	
           std::reverse(R.begin(), R.end());	
           std::reverse(alt.begin(), alt.end());
        std::cout << alt[0] << std::endl;	
        std::cout << alt[num-1] << std::endl;	
// 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, &num, &c__1, nr, &alt[0], zr, 
			&R[0], &q[j][0]  );
	  }
	  else
	  {
		R = ::gLevel1Data.msisO();
// O reduce by 50%
//    std::transform(R.begin(), R.end(), R.begin(),
//        std::bind2nd( std::multiplies<float>(), 0.5   ));
		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);
*/
}