#include "DatabaseHandles.h"
#include "FORTRANarray.hpp"
#include "F77Radtran.h"
#include "LoadLinesFromDB.h"
#include "fortranfunctions.h"
#include "GATS_Utilities.hpp"

#include <iostream>
#include <cassert>

using namespace gatsDBpp;
using namespace GATS_Utilities;

extern "C" void sofie_lbl_model (
	int* lfilt,	int*  lhvy,	int* lout,	int* dbgin,
	int* ibands,	int* nbands,	int* isc,
	int* igas,  	int* mgas,	int* ngas,
	int* mcel,	int* ncel,	int* isl,	int* iel,
	int* mpout,	float*  pout,	float ***rlt,	float ***rut,
	float* re,	float* tearth,	float* albedo,
	float* zt,	float* pt,	float* tt,
	float **qmix,	float ***qg,	
	double  **tau	, float* xlat, float* za

)
{
	const static int crssSize[]  =  {
	      2000000,   1600000,      // s,w channel 1  Ozone          s=strong  w=weak
		7000,   7000,      // s,w channel 2  UV PMC
	       15000,  85200000,	   // w,s channel 3  H2O
	       25000,  25000,      // s,w channel 4  2.7 micron CO2
               10000,  10000,      // s,w channel 5  IR PMC
               60000,  60000,      // s,w channel 6  CH4
	       50500000, 70000,     // s,w channel 7  4.3 micron CO2
	       40000,  40000, 	   // w,s channel 8  NO
	};


	static float dvel = 0.0;	// Doppler shift... 

	static int oldisc = -999;
	static std::vector<int> oldIbands, oldGasIds;
	static int nfpmx = 2840;		// max number points per filter
	static int icalc = 1;  		// transmission only
	static int ic1 = -1;		//beginning cell for integration
	static int ic2 = 1;		//ending cell for integration
    static int iopt = 2;
    static int imr = 2;   // multiple rays
    static int nset=1;
    static int c_one=1;
    static int initIntegrate=1;

	static int mcout, mray, nray;

	static std::string lineTableName = "hitran2004" ;
    static std::string connectionName = "LineData";

	static double vs1;
    static double vs2;

	static FORTRANarray<float> pmass, tbs, tedge, gasb,rl,ru,tb, pb;
	static std::vector<float> filt, coutA;
	static std::vector< std::vector<float> >  Vcntrl;
	static std::vector< std::vector<double> >  Vcrss;


    double avrad, avtran, filtti;


    if(pmass.size() == 0) {
        mray= *mcel;
        pmass = FORTRANarray<float>( *mgas, *mcel, mray, 2);    // static
        tbs = FORTRANarray<float>( *mcel, mray, 2, *mgas);  //static
        tedge=FORTRANarray<float>( *mcel, mray, 2, *mgas);  //static
        gasb = FORTRANarray<float>(*mgas, *mcel);       //static
        ru = FORTRANarray<float>( *mcel, mray, *mgas);
        rl = ru;
        tb = FORTRANarray<float> (*mcel);
        pb = FORTRANarray<float> (*mcel);
	}
//
// determine if the band ids and the atmosphere ID are the same as the previous call
//   if they are different we need to calculate cross sections
//   if they are the same, we can skip to the integrator only



	if( oldisc != *isc  || 
        *nbands != (int)oldIbands.size() || 
	*ngas   != (int)oldGasIds.size() ||
	! std::equal(oldIbands.begin(), oldIbands.end(), ibands) ||
	! std::equal(oldGasIds.begin(), oldGasIds.end(), igas)   ) {
//std::cout << " only once " << std::endl;

/*
std::cout << *xlat << "   " << *re << std::endl;
std::cout << "&&&&&&" << std::endl;
        for(int i= *ncel-1; i>=0; --i )
std::cout << zt[i] << "  " << pt[i] << std::endl;
std::cout << "&&&&&&" << std::endl;
        for(int i= *ncel-1; i>=0; --i )
std::cout << zt[i] << "  " << tt[i] << std::endl;
std::cout << "&&&&&&" << std::endl;

for (int j=0; j< *ngas; j++ ) {

std::cout << igas[j] << std::endl;

for(int i= *ncel-1; i>=0; --i )
 std::cout << zt[i]  << " " <<  qmix[j][i]*1.e6  << std::endl;
std::cout << "&&&&&&" << std::endl;
}
*/


		nray =  *ncel;			//static
		oldisc = *isc;			//static
		oldIbands = std::vector<int>(ibands, ibands+ *nbands);  //static
		oldGasIds = std::vector<int>(igas, igas+  *ngas);  //static
		int mfilt = *nbands * (nfpmx+10);
		filt = std::vector<float>(mfilt);		// static
		Vcntrl.clear();
		Vcrss.clear();

    	int it2 = *ncel *  *ncel;
    	int it1 = it2 + *ncel;
   		mcout = *ncel * 11 + it1 * 6 + it2 * 3;  // static
		double resfilt;

        std::vector<float> amass(pmass.size(), 0.0) ;
	coutA = std::vector<float>(mcout);     // static

// mass apportioning
		::lp_celltb__( mcel, &mray, mgas, 
			mpout, &mcout, ncel, &nray,  
			&nset, &imr, &c_one, ncel, ngas, 
			igas, re, pt, tt, 
            &rlt[0][0][0], &rut[0][0][0], &qmix[0][0],  &qg[0][0][0],  
  		    tb.getArray(), tbs.getArray(), 
			tedge.getArray(), rl.getArray(), 
			ru.getArray(), pb.getArray(), gasb.getArray(), 
			pmass.getArray(),  pout, &coutA[0]   );

// Read in the Filters for the requested bands
		float solar=0.0;		// No solar weighting applied in Lp_filter
		::lp_filter__( lfilt, nbands, &mfilt, 
			nbands, &nfpmx, &solar, &vs1, 
			&vs2, ibands, &resfilt,  &filt[0]);

		DatabaseHandles* handles = GetDatabaseHandles();
		assert(handles->HandleExists(connectionName));
		GATS_DB  *dbConnection = handles->Get(connectionName).get();
// 
		int mbnd = 5830;		// Max number of points for band model
		int mcof = 6;		// Max number of coefficients for each band model
		int mset = 1;          // only one cross section set per gas;

		int mcntrl = ( *ngas << 1) + 6 + *ncel * ( *ngas * 13 + 5 +  *ngas * mray);	
		int nhvy=4+ *ngas * (2+ 11 * mset+ mcof * mbnd * mset);
		std::vector<float> hvy(nhvy);
        Vcntrl.push_back(std::vector<float>(mcntrl)); //, cntrl(mcntrl);
		std::vector<int> gasId(igas,igas+ *ngas);

		for(int ifilt= 1; ifilt <= *nbands; ++ifilt ) {

			::filterlimits_( &vs1, &vs2, &ibands[ifilt-1], nbands, &filt[0] );

			::lp_hvyset__( lout, lhvy, dbgin,
	 			&mcof, &mbnd, &mset, ngas, igas, 
	 			&vs1, &vs2,  &hvy[0]  );

			float wnmax = 50.0; // need to change this
			int mcmix = 1;			// max number of parameters for line mixing
			int mlmix = 5;         // Max number of lines allowed to have mixing
			std::vector<float> afgl;
//std::cerr << "here " << std::endl;
			int nl = LoadLinesFromDB( gasId, vs1, vs2, wnmax, 
				mcmix, mlmix,  afgl,  dbConnection ,lineTableName);
//std::cerr << "here " << std::endl;
			int segment=0;	// don't segment here
			double resmax = 0.;   // Linepak picks its own resolution
			float odmin = -1.0;  // if negative Linepak will use its default value 
			float resfc = 0.3;  	// Resolution factor 
            int initCross=1;
//std::cout << "here I am " << nl <<  std::endl;
			int msiz = crssSize[ ibands[ifilt-1]-1 ];
			Vcrss.push_back(std::vector<double>(msiz) ); // crss(msiz);
			int  done;

			::lp_crossm__(  
				lout,  &iopt, &icalc,
	 			mcel, &mray,  mgas, &mcntrl, &msiz, 
	 			&nl,  &mcmix,  ncel, &nray, 
				ngas, igas, &vs1,  &vs2, 
				&wnmax, &resfc, &resmax, &odmin, 
				&c_one, ncel,    // rays 1 to ncel
				pb.getArray(), tb.getArray(), gasb.getArray(), 
				&amass[0], &segment, &((Vcntrl[ifilt-1])[0]), &afgl[0], &hvy[0], 
				&((Vcrss[ifilt-1])[0]), &initCross, &done, &nset, &imr, 
				mpout,  pout, rl.getArray(), ru.getArray()  );

			afgl.clear();
//initIntegrate=1;

/*
*isl=1;
*iel= *ncel;
*/


			for(int ir = *isl; ir <= *iel; ++ir) {
				::lp_cingrate__( lout, &icalc, &iopt, mcel, 
					&mray,  mgas, ncel, &nray, &initIntegrate, 
					&ir, &ic1, &ic2,
					ngas, igas, &vs1, &vs2, albedo, 
					tearth, tbs.getArray(), 
					tedge.getArray(), pmass.getArray(), pout,
					&((Vcntrl[ifilt-1])[0]), &((Vcrss[ifilt-1])[0]) )	;
				
				::lp_cfilter__( &icalc, nbands,
					&ifilt,  nbands, &nfpmx, &dvel, &avtran, 
					&avrad, &filtti, &((Vcntrl[ifilt-1])[0]), &((Vcrss[ifilt-1])[0]), &filt[0] );

				tau[ifilt-1][ir-1] =  avtran / filtti;
//std::cout << "inside top " << za[ir-1] << " " << qmix[0][ir-1] << " " << ConvertToStringPrec( tau[ifilt-1][ir-1] )  << std::endl;

			}
/*
for(int ir= *isl-1; ir <= *iel-1; ++ir)
    std::cout << za[ir] << "   " << ConvertToStringPrec(tau[ifilt-1][ir])  << std::endl;
    exit(23);
*/

		} // end loop over bands

// subsquent call for a channel atmosphere with an updated Qmix array
// this is the integrator section

	} else {
        ::lp_celltb__( mcel, &mray, mgas,
            mpout, &mcout, ncel, &nray,
            &nset, &imr, isl, iel, ngas,
            igas, re, pt, tt,
            &rlt[0][0][0], &rut[0][0][0], &qmix[0][0],  &qg[0][0][0],
            tb.getArray(), tbs.getArray(),
            tedge.getArray(), rl.getArray(),
            ru.getArray(), pb.getArray(), gasb.getArray(),
            pmass.getArray(),  pout, &coutA[0]   );
/*
		::lp_upmass__( mcel, &mray, mgas,  mpout, ncel, 
			&nray, isl,  iel,  ngas,
			igas,  &qmix[0][0], &qg[0][0][0],  
			gasb.getArray(), pmass.getArray(), pout) ;
*/
//initIntegrate=1;
			
		for(int ifilt= 1; ifilt <= *nbands; ++ifilt ) {   // loop over the bands

            for(int ir = *isl; ir <= *iel; ++ir) {
                ::lp_cingrate__( lout, &icalc, &iopt, mcel,
                    &mray,  mgas, ncel, &nray, &initIntegrate,
                    &ir, &ic1, &ic2,
                    ngas, igas, &vs1, &vs2, albedo,
                    tearth, tbs.getArray(),
                    tedge.getArray(), pmass.getArray(), pout,
                    &((Vcntrl[ifilt-1])[0]), &((Vcrss[ifilt-1])[0]) )   ;

                ::lp_cfilter__( &icalc, nbands,
                    &ifilt,  nbands, &nfpmx, &dvel, &avtran,
                    &avrad, &filtti, &((Vcntrl[ifilt-1])[0]), &((Vcrss[ifilt-1])[0]), &filt[0] );

                tau[ifilt-1][ir-1] =  avtran / filtti ;

//std::cout << "Later On " << za[ir-1] << " " << qmix[0][ir-1] << " " << ConvertToStringPrec( tau[ifilt-1][ir-1] )  << std::endl;
//std::cout << ir << " Later On " << ConvertToStringPrec(filtti) << std::endl;
            }

		}  // loop over bands
	} // end else type
}