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

#include <fstream>
#include <iostream>
#include <cassert>
#include <cmath>
#include <functional>
#include <cstdlib>
/*
extern "C" void calc_sim_exo ( int* nblk,   int*  iblk,
        double* exolock ,  double *simexo, int* fovflag, double* fov_offsets);
*/

using namespace GATS_Utilities;

extern "C" void get_level1_data_replace( int* isignals,  int* nsignals, 
		int* ncel,	float* za, int* pmcflag,	
		double **meas_signal, double* appang, double* lock,
		int* iblk,	int* nch,
		int* ref_band, int* extrap_flag,  int* iexo,
		char*  filename ) 
{

        srand(23000);

        static const float Ratio[10][16] = {
// band   1     2    3   4    5  6    7       8       9    10      11    12        13         14       15     16
        { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.}, //1
//        { 1.719, 0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.},  //2
        { 1.55, 0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,     0.,      0.,       0.,      0.},  //2
	    { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.}, //3
        { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.}, //4
        { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.},  //5
        { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.},  //6
        { 0.,   0.,  0., 0., 0., 0., 0.,     0.,     0.,   0.,  0.,      0.,	  0.,      0.,       0.,      0.},  //7 						
        { 0.,   0.,  0., 0., 0., 0., 0.045 , 0.,     0.,   0.,  0.0835,  0.,      0.0585,  0.0486,   0.,      0.},   // band 8
        { 0.28, 0.,  0., 0., 0., 0., 0.0155, 0.3857, 0.,   0.,  0.0322,  0.0123,  0.025,   0.0187,   0.0085,  0.},   //band 9
        { 0.,   0.,  0., 0., 0., 0., 0.0332, 0.8261, 0.,   0.,  0.0695,  0.0266,  0.0484,  0.0402,   0.,      0.} }; // band 10


	const static int bandsToCorrect[] = { 1, 7, 11, 13 } ;  // sofie bands that get an aerosol correction
								// this is ozone 2.8 CO2, CH4, and 4.3 CO2


    static int c__1 = 1;

    if(isignals[0] == 500) return;   // this is the refraction angles
    std::ifstream inputFile;   //the input file

    std::vector<double> altgrid =  ::gLevel1Data.GetImpactAlts() ;
    assert( ! altgrid.empty() );

//    std::transform( altgrid.begin(), altgrid.end(), altgrid.begin(),
//                    std::bind2nd( std::plus<double>(),  0.5 ) );
//    std::cout << "added alts " << std::endl;


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

//std::cout << " hello " << std::endl;
    inputFile.open(filename, std::ios_base::in);

    if(inputFile.is_open()    ) {   /// read the file
        std::vector<double> A,B;
        double T,Z;
        int nlev,iband;


        inputFile >> nlev ;
        inputFile >> iband;

         double ExoC = ::gLevel1Data.GetExoCounts( iband )  ;

        for(int rr =0 ;  rr < nlev; rr++) {
                inputFile >> T;
                inputFile >> Z;
// if(Z < 97. && Z > 93.) std::cout << ConvertToStringPrec(T) << "  "
//   <<  ConvertToStringPrec(Z) << std::endl;
               A.push_back(T*ExoC ) ;   B.push_back(Z);
//            A.push_back(T  + 2.e-5*( (rand() % 5)-2)  )  ; B.push_back(Z);
        }
        inputFile.close();
//exit(23);
        assert(monotonic_if( B.begin(), B.end(), std::greater<float>()    )  );

         std::cout << " read from signal file " << std::endl;
         std::vector<double> S(num1) ;

        ::linerd_(&B[0], &A[0], &altgrid[0], &S[0], &nlev, &num1 );

        ::gLevel1Data.SetSignal( iband,S ) ;
   }

    std::vector<double> offsetAng(num1), offsetTran(num1);

    std::vector<double> lockV = ::gLevel1Data.LockDownElev();
    std::vector<double> appangV = ::gLevel1Data.GetViewAngle();

    std::vector<double> zaDouble(za,za+ *ncel), tran;

    assert ( altgrid.size() == lockV.size() );

// FOR TESTING
    lockV = std::vector<double>(num1,16.*M_PI/(60.*180.) ); 
//

    ::linerd_(&altgrid[0], &lockV[0], &zaDouble[0], lock, &num1, ncel );
	
    assert ( altgrid.size() == appangV.size() );
    ::linerd_(&altgrid[0], &appangV[0], &zaDouble[0], appang, &num1, ncel );

//    inputFile.open(filename, std::ios_base::in);  // open the file

    for(int i=0; i<  *nsignals; ++i) {
        int isig = isignals[i] ;	

	    if( isig >=1 && isig <=24) {

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

// Apply FOV correction to signals
           
             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 );

 // add a top point to signals so convolution doesn't screw up so much
/*
            double val, ztop = altgrid.front() + 30.; 
              ::linerd_(&altgrid[0], &appangV[0], &ztop, &val, &num1, &c__1 );
             appangV.insert(appangV.begin(), val);

              ::linerd_(&altgrid[0], &tran[0], &ztop, &val, &num1, &c__1 );
              tran.insert(tran.begin(), val);

              ::linerd_(&altgrid[0], &lockV[0], &ztop, &val, &num1, &c__1 );
              lockV.insert(lockV.begin(), val);

              altgrid.insert(altgrid.begin(), val);
              num1++;
*/

              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] );


 //                for(int i = 0; i< num1; ++i) {
 //                  std::cout << ConvertToStringPrec( 1.0- tran[i] ) << "  "  << altgrid[i] << std::endl;
 //                } 

  //            std::cout << "&"  << std::endl;
  //             std::cout << "# convolved"  << std::endl;

  //               for(int i = 0; i< num1; ++i) {
  //             std::cout << ConvertToStringPrec( 1.0-convtau[i] ) << "  "  << altgrid[i] << std::endl;
 //                 }


 //             std::cout << "&"  << std::endl;
 //              std::cout << "# diff"  << std::endl;

        std::vector<double> diffs, tauRaw;
         for(int i= 0; i < num1; ++i) {
           diffs.push_back(  tran[i]/ convtau[i]     ) ;
//           std::cout << ConvertToStringPrec( diffs.back()  ) << "  "  << altgrid[i] << std::endl;
         }



  //            std::cout << "&"  << std::endl;
  //             std::cout << "# the Raw"  << std::endl;

         for(int i= 0; i < num1; ++i) {
           tauRaw.push_back( tran[i]*diffs[i]         ) ;
//std::cout << ConvertToStringPrec( 1.0- tauRaw.back()   ) << "  "  << altgrid[i] << std::endl;

         }

 ::gLevel1Data.SetSignal( isig,tauRaw ) ;
 tran = tauRaw;
                 

//                   exit(23);
        
// end

	    } else if(isig > 16 ) {
       		tran = ::gLevel1Data.GetSignal( isig );
   	    	int b =  2*(isig-17)+1;
   		    if(b == 5 || b == 15) { b++; }   // for H2O and NO channels strong band is an even number
   		    if( *iexo ) {
//std::cout << " this is the strong " << b << std::endl; exit (23) ;
     		     std::transform( tran.begin(),tran.end(), tran.begin(),
        		std::bind2nd( std::divides<double>(),  ::gLevel1Data.GetExoCounts( b )  ) );

  		     } else { 
     		  std::vector<double> hh = ::gLevel1Data.GetSignal( b );
    		      std::transform( tran.begin(),tran.end(), hh.begin(), tran.begin(), std::divides<float>() );
   		     }
	     } 
	} else if(isig == 25)  {
   	     tran = ::gLevel1Data.SunSensorTransmission();
	} else {
// unknown signal id
	     std::cerr << "Unknown signal ID  " << isig << std::endl;
	     exit(23);	
	}

	assert(tran.size() == altgrid.size() );
	int* j = std::find(iblk, iblk+ *nch, isig);

	if( j != iblk+ *nch) {

// apply aerosol correction to V signals
	   std::map<int, h2>::const_iterator mapIter = ::gRefAero.find( *ref_band ) ;
 
	   if(*pmcflag && *extrap_flag && mapIter != ::gRefAero.end() &&
	   *std::find( bandsToCorrect, bandsToCorrect+sizeof(bandsToCorrect)/sizeof(bandsToCorrect[0]), isig ) == isig) {
                float rat = Ratio[*ref_band-1][isig-1];
		if( *ref_band == 2 && gLevel1Data.EventNo() < 27886 &&  lockV.front() > 14.*M_PI/10800. ) rat = 0.0;
std::cout << " ratio is " << rat << std::endl;
		if(rat <=0.0 ) goto SKIP;
		
                
		std::vector<double> zz = (mapIter->second).zimpact;
		std::vector<double> rr = (mapIter->second).transmittance;

	        double val;
	        int num2=(int)rr.size();

//std::cout << " here is the num size " << num2 << std::endl;

	        std::transform(rr.begin(),rr.end(), rr.begin(), 
		    std::bind2nd(std::ptr_fun( pow), rat ) );
                
//		for(int mm=0; mm<   (int)altgrid.size(); ++mm)
//		  std::cout << ConvertToStringPrec(tran[mm]) << "  " << altgrid[mm] << std::endl;
//		 std::cout << "&" << std::endl;
		 
                for(int i = 0; i< (int)altgrid.size(); ++i) {
                     if(altgrid[i] <= zz.front() && altgrid[i] >= zz.back() )  {
                         ::linerd_(&zz[0], &rr[0], &altgrid[i], &val, &num2, &c__1 );
                           tran[i] /= val;
                      }
                 }
                
//		for(int mm=0; mm<   (int)altgrid.size(); ++mm)
//		  std::cout << ConvertToStringPrec(tran[mm]) << "  " << altgrid[mm] << std::endl;
//		 std::cout << "&" << std::endl;
//		  exit(23);				
	   }   // apply correction

        }  // if

SKIP:


     ::linerd_(&altgrid[0], &tran[0], &zaDouble[0], &meas_signal[i][0],
        &num1, ncel );


    }   //channel loop

}