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

#include <iostream>
#include <numeric>
#include <cmath>

Level1Data gLevel1Data;

using namespace GATS_Utilities;

extern "C" void read_l1_event( int* eventNo, float *re, float *elat, 
	float *elon, char* datestr, int* nl_hi, int* iorbit, int* status,
	double* exolock, char* l1_ver, float* zoffsun, float* rse,
	double** FOVOffsets ) 
{
	static int c__1 = 1;

	if(::gLevel1Data.Read( *eventNo, l1_ver, FOVOffsets ) ) {
		*re = (float)::gLevel1Data.RadiusCurvature();
		*elat = ::gLevel1Data.Lat83();

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

       	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);
*/
			*exolock  = lockV.front() ;

//           std::cout << lockV.front() << "  " << *exolock << std::endl;  exit(23);

// 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>() );
		std::vector<double>::iterator 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;
		unsigned i1 = std::distance(solarExt.begin(), iter1);
		*zoffsun = (float)altgrid[i1] ;

		*rse = gLevel1Data.EarthSunDistance() ;

//std::cout << " read earth " << *rse <<  " off sun " << *zoffsun << std::endl;  exit(23);

//exit(23); 
                *status = 0;
		std::string t = ::gLevel1Data.StartTime()
			+ " (" + ConvertToStringPrec( *elat,4) 
      			+ ", "+ ConvertToStringPrec( *elon,5) + ")";
		strncpy(datestr, t.c_str(), t.size() );
		std::cout << "**  read data for Event No " 
			<< ::gLevel1Data.EventNo() << " " << t << std::endl;
//*status = 1;

// do FOV Removal on SIgnals
   std::cout << "Doing 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(  tran[i]/ convtau[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( isig )  ) )  ;
   ::gLevel1Data.SetSignal( isig,tauRaw ) ;

}

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(  tran[i]/ convtau[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 ) ;



        } else {
		*status = 1;
	}
}