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

#include <algorithm>
#include <iostream>
#include<cmath>

using namespace GATS_Utilities;

extern "C" void make_dv_extinction (
	int* iblk,	int* nch,
	int* idg,  	int* ng,	
	int* nl, 	int* ml,
	int **idg_b,	int **irad_b,	int **icor_b,
	int* itab,	int* mgs,	
	int* ntb,	int *ntt,	int *nst,
	float*  pout,	float ***rlt,	float ***rut,
	float* re,	float* tearth,	float* albedo,
	float **qmix,	float ***qg,	
	float ***tbl,	float *tout,	float **sout,
	double ***emr,  double ***rmr,  double ***esr,
	double ***rsr,	float ***bsr, float* zt,
	float* zapr,   float** qapr, float*** qgapr,
    int* igdo,   int* ngdo, int* ichando, int* nchn )
	
{

        static int c__1= 1;
	FORTRANarray<float> qmixB, qgB, rltB, rutB;

	 qmixB = FORTRANarray<float>(*ml, *mgs);
     qgB = FORTRANarray<float>(2, *ml, *mgs);
	 rltB =FORTRANarray<float>(*ml, *mgs);
	 rutB = FORTRANarray<float>(*ml, *mgs);
    
   FORTRANarray<float> ext(2, *ml);
   FORTRANarray<double> ttran(2,*ml);

//*nch = 1;

	for(int ich = 0; ich< *nch; ++ich) {   // loop over the bands
		int idch = iblk[ich]-1;
//std::cout << *isl << " " << *iel << std::endl;
		int ngs = 0;
		std::vector<int>  iradB, icorB;	
		for(int ig = 0; ig< *mgs; ++ig) {
			if(idg_b[idch][ig] == 0 )  break;
//std::cout <<" idch " <<  idch << std::endl;
			++ngs;
			int* j = std::find(idg, idg+ *ng, idg_b[idch][ig] );
			int ik = (int)std::distance(idg, j);
			if( ik < *ng) {
//std::cout << "# pre " << ik << "  " << ig << "  "  << *ng << std::endl;
				iradB.push_back( irad_b[idch][ig] );
				icorB.push_back( icor_b[idch][ig] );
				std::copy( &qmix[ik][0], &qmix[ik][*nl], qmixB.getArray(1,ig+1) );
//std::fill(qmixB.getArray(1,ig+1), qmixB.getArray(*nl,ig+1) +1, 1.e-4);  
				std::copy( &rlt[ich][ik][0], &rlt[ich][ik][*nl], rltB.getArray(1,ig+1) );
                std::copy( &rut[ich][ik][0], &rut[ich][ik][*nl], rutB.getArray(1,ig+1) );
				std::copy( &qg[ik][0][0], &qg[ik][0][0]+ 2* *nl, qgB.getArray(1,1,ig+1) ); 

			} else {
				iradB.push_back( 0 );
				icorB.push_back( 0 );
			}
		}
		if( ngs == 0) continue;

//std::cout <<"# band " << iblk[ich] << std::endl;

for(int ir= 1; ir <= *nl ; ++ir) {

		int ic1= -ir, ic2= ir; // ic2= ir;
//   std::cout << "qmix " << qmixB(ir, 1) << " " << rltB(ir,1) << "  " << rutB(ir,1) << 
//   "  " << iradB[0] << " " << icorB[0] << std::endl;
   
		::mp_mega__( ntb,  nst,  ntt,  mgs,    ml,    nl,   
			&ir, &ir,   &ic1,  &ic2,   itab,  &ngs,  &idg_b[idch][0],
			qmixB.getArray(), qgB.getArray(), 
			rltB.getArray(), rutB.getArray(),
			albedo, tearth,   &iradB[0],  &icorB[0],   &emr[0][0][0],
			&rmr[0][0][0],  &esr[0][0][0],  &rsr[0][0][0],
			&bsr[0][0][0],  &tbl[ich][0][0],   pout,
			&sout[0][0], tout  );

		double V  = (double)1.0 - emr[0][0][ir-1] ;

//std::cout << zapr[ir-1] << "   " << V << std::endl;
		ext.setVal( -std::log(V)/ ::pathlength_(&ir,&ic1,&ic2,pout), ich+1, ir ) ;
//std::fill(&emr[0][0][0],  &emr[0][0][ir-1]+ *ml * *mgs, 0.0);
/*
        ::mp_mega__( ntb,  nst,  ntt,  mgs,    ml,    nl,
            &ir, &ir,   &c__m1,  &c__1,   itab,  &ngs,  &idg_b[idch][0],
            qmixB.getArray(), qgB.getArray(),
            rltB.getArray(), rutB.getArray(),
            albedo, tearth,   &iradB[0],  &icorB[0],   &emr[0][0][0],
            &rmr[0][0][0],  &esr[0][0][0],  &rsr[0][0][0],
            &bsr[0][0][0],  &tbl[ich][0][0],   pout,
            &sout[0][0], tout  );

        ttran.setVal(   emr[0][0][ir-1],ich+1,ir) ;
 */       
			
//if(*iel - *isl > 10) std::cout << "From Mega  " << za[ir] << " " << qmixB(ir+1,1) << " " << ConvertToStringPrec( tau[ich][ir] )  << std::endl;
}


	} // loop over bands
// now we need to create the difference extinction profile and reset the retrieval values

	for(int ir= 1; ir <= *nl ; ++ir) {
	
		qmix[0][ir-1] = ext(1,ir) - ext(2,ir);   //band 3 - band 4 difference extinction (interferents)

		qmix[1][ir-1] = 1.e-11;  // aerosol   difference extinction


	}

std::fill (&qgapr[0][0][0], &qgapr[0][0][0]+ 2* *ml * *mgs,  0.0 );
std::fill (&rlt[0][0][0], &rlt[0][0][0]+ *ml * *mgs,  1.0 );
std::fill (&rut[0][0][0], &rut[0][0][0]+ *ml * *mgs,  1.0 );

::kp_intprof_noxtrp__( nl, nl, &c__1,
      nl,  zt,  zapr, &qmix[0][0], &qapr[0][0] );
::kp_intprof_noxtrp__( nl, nl, &c__1,
      nl,  zt,  zapr, &qmix[1][0], &qapr[1][0] );

    int idch = iblk[0]-1;
//std::cout << iblk[0] << std::endl; exit(23);
    *nchn = 1;
    ichando[0] = iblk[0]; 

	*ng = 2;  // this was changed
	*nch = 1;
	idg[0] = idg_b[idch][0] = 9001;// 9001;
	idg[1] = idg_b[idch][1] =9000;
	
	irad_b[idch][0] = 4;
	irad_b[idch][1] = 4;
irad_b[idch][2] = 0;
	
	icor_b[idch][0] = 3;
	icor_b[idch][1] = 2;
icor_b[idch][2] = 0;
	

    *ngdo = 1;
    igdo[0] = 9000;

/*
std::cout << "# band " << iblk[0] << std::endl;
for(int t = 0; t< *nl; ++t) 
std::cout << ConvertToStringPrec( ttran(1,t+1)) <<  "   " << zapr[t] <<   "   " << qapr[0][t] << std::endl;
exit(23);

std::cout << "&" << std::endl;
*/

/*
for(int t = 0; t< *nl; ++t)
std::cout << ConvertToStringPrec( ttran(2,t+1)) <<  "   " << zapr[t] <<  std::endl;
exit(23);
*/

//std::cout << "# now simulations of extinction only" << std::endl;

}