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

#include <numeric>
#include <cmath>
#include <cassert>
#include <iostream>
#include <ctime>

#include "Level1Data.h"

float Band7Ratio ( const float& band9, const float& band10 )
{

	const static float enoi = 1.e-7;

	return (band9 > enoi && band10 > enoi ) ? 0.0206571-0.00238340*(band9/band10) : 0.0 ;
}

using namespace GATS_Utilities;

extern "C" void extrapolate_pmc_extinction(
	float* zapr, 	float	**qapr, 	int* napr, 	int* igs, 
	int* ngg,		int* ichando,	int* nchn,	int* extrap_flag,
	int * refband  )
{
	static const float Ratio[3][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.0401, 0.,     0.,   0.,  0.0835,  0.,      0.0585,  0.0486,   0.,      0.},   // band 8
 	{ 0., 0.,  0., 0., 0., 0., 0.0155, 0.3857, 0.,   0.,  0.0322,  0.0123,  0.0226,  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

	if( *extrap_flag == 0) return;
	assert( (*refband>= 8 && *refband <=10) || ( (*refband== 30 || *refband == 40) && ichando[0] == 7 ) );
/*
	int eventNo = gLevel1Data.EventNo();

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

        int idoy = 1 + timestruct.tm_yday;
// check that it is the pmc season   odd events are rises (North hemisphere)  Sets are even (south hemisphere)
	if( ( eventNo % 2 == 1 && (idoy < 135 || idoy > 258)  )   ||  //rise events 15-may through sept 15
            ( eventNo % 2 == 0 && (idoy > 60  && idoy < 305)  ) )  // set events 01-nov - 01-mar
	return;
*/
	int idg, nr;
	std::map<int, h2>::iterator Ziter, Ziter2, Ziter3;

	if(*refband <= 10 ) {
		Ziter = ::gRefAero.find( *refband);
		if( Ziter == ::gRefAero.end() ) {
			std::cout << " Warning Aerosol Retrieval for Band "<< *refband <<" has not been completed" << std::endl;
			return;
		}
 	} else if(*refband == 30) {
		Ziter = ::gRefAero.find( 9 );     // Ziter is the band 9 extinction
		Ziter2 = ::gRefAero.find( 10 );
		if(Ziter == ::gRefAero.end() || Ziter2 == ::gRefAero.end() ) {
			std::cout <<  "Can not perform the Band 9/10 ratio correction " << std::endl;
			return;
		}
	} else {
		Ziter = ::gRefAero.find( 9 );     // Ziter is the band 9 extinction
		Ziter2 = ::gRefAero.find( 10 );
		Ziter3 = ::gRefAero.find( 8 );
		if(Ziter == ::gRefAero.end() || Ziter2 == ::gRefAero.end() ||  Ziter3 == ::gRefAero.end() ) {
			std::cout <<  "Can not perform the Band 8,9,10 extinction average correction " << std::endl;
			return;
		}
	}

	nr = (int)(Ziter->second).zalt.size();
	
	float *i1 = std::find_if( zapr, zapr+ *napr, std::bind2nd( std::less<float>(),95.0) );
	if( i1 == zapr+ *napr || *i1 < 70.0 ) return;

	float *i2 = std::find_if( zapr, zapr+ *napr, std::bind2nd( std::less<float>(),70.0) )-1;


	int isl = (int)std::distance(zapr, i1)+1;  // FORTRAN indicies
	int iel = (int)std::distance(zapr, i2)+1;  // FORTRAN indecies

	for(int ich = 0; ich< *nchn; ++ich) {
		for(int i = 0; i< *ngg; ++i) {
			if( igs[i] > 10000) {
				idg = igs[i] / 10000;	
			} else {
				idg = igs[i] ;
			}
			if( idg >= 9000) {
				std::vector<float> RefBand( *napr), Band10(*napr, 0.0) ;

     				::kp_intprof_noxtrp__( napr, &nr, &isl, &iel,    // for refband =30 or 40 this is band 9
						&( (Ziter->second).zalt[0]),
			      			zapr , 
               		 			&( (Ziter->second).extinc[0]) ,     
						&RefBand[0]  );

				if( *refband == 30 || *refband == 40 ) {
					int npts = (int)(Ziter2->second).zalt.size() ;
      					::kp_intprof_noxtrp__( napr, &npts, &isl, &iel,    // this is band 10 
						&( (Ziter2->second).zalt[0]),
			      			zapr, 
               		 			&( (Ziter2->second).extinc[0]) ,     
						&Band10[0]  );
					if (*refband == 40) {
						std::transform(RefBand.begin(),RefBand.end(), RefBand.begin(),     // Band 9 extrapolated to band 7
							std::bind2nd( std::multiplies<float>(), Ratio[1][6]   ) );

						std::transform(Band10.begin(),Band10.end(), Band10.begin(),     // Band 10 extrapolated to band 7
							std::bind2nd( std::multiplies<float>(), Ratio[2][6]   ) );

						npts = (int)(Ziter3->second).zalt.size() ;
						std::vector<float> Band8(*napr, 0.0);
						::kp_intprof_noxtrp__( napr, &npts, &isl, &iel,    // this is Band 8
							&( (Ziter3->second).zalt[0]),
			      				zapr, 
               		 				&( (Ziter3->second).extinc[0]) ,     
							&Band8[0]  );

						std::transform(Band8.begin(),Band8.end(), Band8.begin(),     // Band 8 extrapolated to band 7
							std::bind2nd( std::multiplies<float>(), Ratio[0][6]   ) );
/*
std::cout << "#Band 8 correction" << std::endl;
for(int il =isl-1; il < iel; ++il ) {
	std::cout <<Band8[il]  << " " << zapr[il] << std::endl;
}
std::cout << "&" << std::endl;

std::cout << "#Band 9 correction" << std::endl;
for(int il =isl-1; il < iel; ++il ) {
        std::cout <<RefBand [il]  << " " << zapr[il] << std::endl;
}
std::cout << "&" << std::endl;

std::cout << "#Band 10 correction" << std::endl;
for(int il =isl-1; il < iel; ++il ) {
        std::cout <<Band10[il]  << " " << zapr[il] << std::endl;
}
std::cout << "&" << std::endl;
*/
						std::transform(RefBand.begin(),RefBand.end(), Band8.begin(), RefBand.begin(),
							std::plus<float>() );
						std::transform(RefBand.begin(),RefBand.end(), Band10.begin(), RefBand.begin(),
							std::plus<float>() );
						std::transform(RefBand.begin(),RefBand.end(), RefBand.begin(), // RefBand is now average extinction
							std::bind2nd( std::divides<float>(),3.0) );   // from extrapolated bands 8,9, and 10
/*
std::cout << "#Mean correction" << std::endl;
for(int il =isl-1; il < iel; ++il ) {
        std::cout <<RefBand[il]  << " " << zapr[il] << std::endl;
}
std::cout << "&" << std::endl;
*/

					}
				}	

				for(int il = isl-1; il < iel; ++il ) {
					float factor;
					if( *refband == 40 ) {
				 		factor = 1.0;
					} else if( *refband == 30) {
						factor = Band7Ratio( RefBand[il],Band10[il]) ;
					} else {
						factor = Ratio[ *refband-8][ichando[ich]-1] ;
					}
					if(factor < 1.e-4) continue;
//		std::cout <<	zapr[il] << " " << factor << " " << RefBand[il] <<  std::endl;	
					qapr[i][il] = factor*RefBand[il];
				}
			}
		}
	}
}