#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_extinction2 (
	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,
	float* za, double** meas_sig )
	
{
//std::cout << *nch << std::endl;

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

	for(int ich = 0; ich< *nch; ++ich) {   // loop over the bands
		int idch = iblk[ich]-1;
		int ngs = 0;
		std::vector<int>  iradB, icorB;	
		for(int ig = 0; ig< *mgs; ++ig) {
			if(idg_b[idch][ig] == 0 )  break;
			++ngs;
			int* j = std::find(idg, idg+ *ng, idg_b[idch][ig] );
			int ik = (int)std::distance(idg, j);
			if( ik < *ng) {
				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;
/*
for(int kk=0; kk< ngs; ++kk)
   std::cout << idch << "  " << iradB[kk] << "  " << icorB[kk] << " " << idg_b[idch][kk] << std::endl;

for(int mm=1; mm<= *nl; ++mm ) {
 
    for(int jj=1; jj<= ngs; ++jj) 
      {
         std::cout << qmixB(mm, jj) << "  " ;
      }
  std::cout << zt[mm-1] << std::endl;
}

std::cout << " next turn " << std::endl;
*/
        for(int ir= 1; ir <= *nl ; ++ir) {
	    	int ic1= -1, ic2= 1; 
			::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] ;
			ext.setVal( V, ich+1, ir ) ;
		}  // loop over layers

	} // loop over bands

// now we need to create the difference extinction profile and reset the retrieval values


//for(int ir= 1; ir <= *nl ; ++ir) {
//std::cout << ConvertToStringPrec(meas_sig[0][ir-1]) << " " << za[ir-1] << std::endl;
//}
//std::cout << "&" << std::endl;

/*
for(int ir= 1; ir <= *nl ; ++ir) {
std::cout << ConvertToStringPrec( ext(1,ir) ) << "  " << ConvertToStringPrec( ext(2,ir)  ) << " " << za[ir-1] << std::endl;
}
std::cout << "&" << std::endl;

*/


	for(int ir= 1; ir <= *nl ; ++ir) {
//		qmix[1][ir-1] = ext(1,ir) - ext(2,ir);   //band 3 - band 4 difference extinction (Rayleigh + gases)
		qmix[0][ir-1] = 1.e-10;  // residual extinction retrieved
// now remove rayleigh from measurement
//      std::cout << meas_sig[0][ir-1] << "  " <<  meas_sig[0][ir-1] - ( ext(2,ir) - ext(1,ir) ) << "  " << za[ir-1] << std::endl;
       meas_sig[0][ir-1] -= ( ext(2,ir) - ext(1,ir) );
	}
//exit(23);

/*     
for(int ir= 1; ir <= *nl ; ++ir) {
std::cout << ConvertToStringPrec(meas_sig[0][ir-1]) << " " << za[ir-1] << std::endl;
}
exit(23);
*/


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


/*
	std::vector<float> zN (&zt[0], &zt[*nl] );
   
    for(int j = 1; j< *nl; ++j ) {
      zN[j] +=  0.25*(  zt[j-1]-zt[j] );
    }
*/

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

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

	*ng = 1;  // this was changed
	*nch = 1;
//idg[1] = idg_b[idch][0] = 0;
	idg[0] = idg_b[idch][0] =9000;
	idg_b[idch][1] = 0;
	
	irad_b[idch][0] = 4;
//irad_b[idch][1] = 4;
	irad_b[idch][1] = 0;
	
	icor_b[idch][0] = 3;
//icor_b[idch][1] = 2;
	icor_b[idch][1] = 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 << qapr[0][t] << "   " << qapr[1][t]  <<  "   " << zapr[t] <<  std::endl;
exit(23);
*/

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

}