#include "Level1Data.h"
#include "GATS_Utilities.hpp"
#include "fortranfunctions.h"
#include "sofiefov.h"

#include <iostream>
#include <numeric>
#include <cmath>
#include <fstream>
#include <cassert>
#include <cstdlib>

//Level1Data gLevel1Data;

using namespace GATS_Utilities;

extern "C" void get_average_l1_data( char* flist, float *re, float *elat, 
	float *elon, char* datestr, int* nl_hi, int* iorbit, int* status,
	double* exolock, char* l1_ver, float* zoffsun, int* pointing65,
	float* rse, int* fovflag, int* poorman, int* ievent ) 
{

static int c__1 = 1;
	std::ifstream Infile(flist);
	assert(Infile);
	
	std::istream_iterator<int> b(Infile),e;
	std::vector<int> events(b,e);
	Infile.close();
	
	assert(! events.empty() );
	std::cout << "**** Averaging Events  *****" << std::endl;
	for(unsigned i = 0; i< events.size() ; ++i) 
	{
		std::cout << events[i] << std::endl;
	}
	
	if(::gLevel1Data.AverageList( events, l1_ver ) ) {
		*re = (float)::gLevel1Data.RadiusCurvature();
		*elat = ::gLevel1Data.Lat83();

		*elon = ::gLevel1Data.Lon83() ;
		*iorbit = ::gLevel1Data.OrbitNo();
		*ievent = ::gLevel1Data.EventNo();

       	std::vector<double> altgrid =  ::gLevel1Data.GetImpactAlts() ;
        std::vector<double> lockV = ::gLevel1Data.LockDownElev();  // radians

 		*nl_hi = (int) altgrid.size();

// Determine the Exo Lockdown...   We're using an average between 120 and 110 km
	       	std::vector<double>::iterator iter1 = std::find_if(altgrid.begin(),altgrid.end(),
                        std::bind2nd(std::less_equal<double>(),120.) );
        	std::vector<double>::iterator iter2 = std::find_if(altgrid.begin(),altgrid.end(),
                        std::bind2nd(std::less_equal<double>(),110.) );
        	unsigned i1 = std::distance(altgrid.begin(), iter1);
        	unsigned i2 = std::distance(altgrid.begin(), iter2);
        	*exolock  = std::accumulate(lockV.begin()+i1, lockV.begin()+i2+1,0.0)/ (double)(i2-i1+1);

// Determine the altitude (zoffsun) when the Lockdown is within 2.5 arcminutes of the bottom edge.
		std::vector<double> solarExt = ::gLevel1Data.SolarExtent();

		std::transform(solarExt.begin(), solarExt.end(), lockV.begin(), solarExt.begin(),
			std::minus<double>() );
		iter1 = std::find_if( solarExt.begin(), solarExt.end(),
			std::bind2nd(std::less_equal<double>(), M_PI* 2.5/ (60.* 180.) ) ); // 2.5 arcmin from bottom edge
		if( iter1 == solarExt.end() )  --iter1;
		i1 = std::distance(solarExt.begin(), iter1);
		*zoffsun = (float)altgrid[i1] ;

		*rse = gLevel1Data.EarthSunDistance() ;

/*
for(int bt = 1; bt< 17; ++bt) {
  std::vector<double> rmt = ::gLevel1Data.GetSignal(bt);


std::cout << "# band " << bt << "  " << ::gLevel1Data.GetExoCounts(bt) <<std::endl;
  for(int i =0; i < altgrid.size() ; ++i)
     std::cout << ConvertToStringPrec(rmt[i])
     << "  "  << altgrid[i] << std::endl;
std::cout << "&" << std::endl;

}

exit(23);

*/


if (! *fovflag && poorman)
 {
// do FOV Removal on SIgnals
   std::cout << "Doing poorman FOV correction " << std::endl;
   std::vector<double> appangV = ::gLevel1Data.GetViewAngle();

   int num1 = (int)appangV.size();

   for(int isig = 1; isig<= 16; ++isig ) {

     std::vector<double> tran = ::gLevel1Data.GetSignal( isig );
     std::transform( tran.begin(),tran.end(), tran.begin(),
        std::bind2nd( std::divides<double>(),  ::gLevel1Data.GetExoCounts( isig )  ) );

     std::vector<double> fovdat,fovpos;

     std::vector<double> elevSldc = ::gLevel1Data.ScanElev();
     int numsol = (int)elevSldc.size();
     std::vector<double> sldc = ::gLevel1Data.ScanData(isig );

     fovpos = ::gFOVdata.get_fov_elev( isig );
     fovdat = ::gFOVdata.get_fov_response(isig );

     int numfov = fovdat.size();

     std::transform(fovpos.begin(),fovpos.end(),fovpos.begin(),
     std::bind2nd( std::multiplies<double>(), M_PI/(60.*180.) )  );

     double simexo;
     ::calc_simexo__(&fovdat[0], &fovpos[0], &numfov, &sldc[0],
                     &elevSldc[0], &numsol,  &lockV[0] , &simexo );

     std::vector<double> convtau(num1);

     ::conv_fov_atm_sun__(
                   &num1,     &c__1,   &num1,
                   &appangV[0],  &appangV[0], &lockV[0],  &tran[0],
                   &fovdat[0],   &fovpos[0],   &numfov,
                   &sldc[0],   &elevSldc[0],    &numsol,
                   &simexo,  &convtau[0] );

     std::vector<double> diffs, tauRaw;
     for(int i= 0; i < num1; ++i) {
           diffs.push_back(   convtau[i] -tran[i]    ) ;
     }

     for(int i= 0; i < num1; ++i) {
           tauRaw.push_back( tran[i] - diffs[i]   ) ;

     }

/*
if(isig >=2 && isig <=4) {
   std::cout << "# band " << isig << std::endl << "@with g0" << std::endl;
 
  for(int i =0; i < num1; ++i)
     std::cout << ConvertToStringPrec(1.0-tran[i])
     << "  "  << altgrid[i] << std::endl;
std::cout << "&" << std::endl;
  for(int i =0; i < num1; ++i)
     std::cout << ConvertToStringPrec(1.0-tauRaw[i])
     << "  "  << altgrid[i] << std::endl;
  std::cout << "&" << std::endl << "@with g1" << std::endl;

  for(int i =0; i < num1; ++i)
     std::cout << ConvertToStringPrec(
       100.*
   (tran[i]-tauRaw[i]) / ( 0.5 * ( 2.0-tauRaw[i]-tran[i] ) ) )
     << "  "  << altgrid[i] << std::endl;
  std::cout << "&" << std::endl;
}
*/



     std::transform(tauRaw.begin(),tauRaw.end(),tauRaw.begin(),
     std::bind2nd( std::multiplies<double>(),  ::gLevel1Data.GetExoCounts( isig )  ) )  ;

     ::gLevel1Data.SetSignal( isig,tauRaw ) ;

}

//exit(23);

int isig = 18;
     std::vector<double>        tran = ::gLevel1Data.GetSignal( isig );
     std::transform( tran.begin(),tran.end(), tran.begin(),
     std::bind2nd( std::divides<double>(),  ::gLevel1Data.GetExoCounts( 3 )  ) );

     std::vector<double> fovdat,fovpos;

     std::vector<double> elevSldc = ::gLevel1Data.ScanElev();
     int numsol = (int)elevSldc.size();
     std::vector<double> sldc = ::gLevel1Data.ScanData(3 );

     fovpos = ::gFOVdata.get_fov_elev( 3 );
     fovdat = ::gFOVdata.get_fov_response(3 );

     int numfov = fovdat.size();

     std::transform(fovpos.begin(),fovpos.end(),fovpos.begin(),
     std::bind2nd( std::multiplies<double>(), M_PI/(60.*180.) )  );

     double simexo;
     ::calc_simexo__(&fovdat[0], &fovpos[0], &numfov, &sldc[0],
                     &elevSldc[0], &numsol,  &lockV[0] , &simexo );

     std::vector<double> convtau(num1);

     ::conv_fov_atm_sun__(
                   &num1,     &c__1,   &num1,
                   &appangV[0],  &appangV[0], &lockV[0],  &tran[0],
                   &fovdat[0],   &fovpos[0],   &numfov,
                   &sldc[0],   &elevSldc[0],    &numsol,
                   &simexo,  &convtau[0] );
     std::vector<double> diffs, tauRaw;
     for(int i= 0; i < num1; ++i) {
           diffs.push_back(   convtau[i]   - tran[i]   ) ;
     }

     for(int i= 0; i < num1; ++i) {
           tauRaw.push_back( tran[i]-diffs[i]   ) ;
     }


     std::transform(tauRaw.begin(),tauRaw.end(),tauRaw.begin(),
     std::bind2nd( std::multiplies<double>(),  ::gLevel1Data.GetExoCounts( 3 )  ) )  ;
     ::gLevel1Data.SetSignal( isig,tauRaw ) ;


 }


/*
std::cout << " read earth " << *rse <<  std::endl;

std::cout << "off sun " << *zoffsun << std::endl;

std::cout << " exo lock " << *exolock << std::endl;

std::cout << " band 2 exo counts " << ::gLevel1Data.GetExoCounts(2) << std::endl;

std::cout << " band 2 Attenuator " << ::gLevel1Data.Attenuator(2) << std::endl;


std::vector<float>  zz = ::gLevel1Data.Alt();
std::vector<float>  pp = ::gLevel1Data.Pres();
std::vector<float>  tt = ::gLevel1Data.Temp();

for (int i=0; i< zz.size() ; ++i ) {

std::cout << zz[i] << "  " << pp[i] << "   " << tt[i] << std::endl;

}
exit(23);



std::vector<double> Sig = ::gLevel1Data.GetSignal(2);   
for(int i=0; i< Sig.size() ; ++i ) {

 std::cout << Sig[i] << "    " << altgrid[i] << std::endl;

}

exit(23); 
*/

 *status = 0;

//*status = 1;

 } else {
        *status = 1;
 }

//std::cout << "finished averages " << std::endl; exit(23);
	
	
}