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

#include <iostream>

using namespace GATS_Utilities;
extern "C" void extend_tran_profile ( int* iblk,  int* nch, 
		int* isl,     int*   iel,    int*  ncel,
		int* layer_start,   int* layer_stop,  int* bottom_up,
		float* za,  double** tau, int* fov_flag, int* ifov_loop  )
{

	static std::vector< std::vector<double> > old_tran;
	double delta;

//   std::cout <<  *isl << "  " << *iel << "  " << *layer_start << "   " << *layer_stop <<  " " << *fov_flag << std::endl; 
//std::cout << " ib " << iblk[0] << std::endl;
//    if(iblk[0] == 13 ) {
//std::cout << " FOV flag is " <<  *fov_flag << std::endl;
  // exit(23);
//}
 
//    if( *ifov_loop > 1) *fov_flag = 1; 
	if( ! *fov_flag) return;  // don't need to do this if no FOV convolution
	
// This will block will execute only on the first call for a channel
//  std::cout << " tau " << tau[0][0] << std::endl; exit(23);
	if( std::find_if( &tau[0][0], &tau[0][0]+ *ncel, 
     std::bind2nd(std::less<double>(), -900. ) ) != &tau[0][0]+ *ncel )  {
		old_tran.clear();
//std::cout << " here " << std::endl; exit(23) ; 


		std::vector<double> altgrid =  ::gLevel1Data.GetImpactAlts() ;
		int num1 = (int)altgrid.size();
		std::vector<double> zaDouble(za,za+ *ncel), tran( *ncel) ;
		for( int ich = 0; ich< *nch; ++ich ) {

                  std::vector<double> meas_tran = ::gLevel1Data.GetSignal( iblk[ich] );
                  std::transform( meas_tran.begin(),meas_tran.end(), meas_tran.begin(),
                    std::bind2nd( std::divides<double>(),  ::gLevel1Data.GetExoCounts(  iblk[ich] )  ) );
			::linerd_(&altgrid[0], &meas_tran[0], &zaDouble[0], &tran[0], &num1, ncel );


//             std::cout << ConvertToStringPrec(1.-tau[ich][i]) << "   " << zaDouble[i] << std::endl;
               
		 
//    std::cout << "&" << std::endl;
//   for(int i= 0; i< *ncel ; ++i) {

//             std::cout << ConvertToStringPrec(1.0-tran[i]) << "   " << zaDouble[i] << std::endl;
//        }
// std::cout << "&" << std::endl;


// Tau is the simulated transmission
// scale the measured above and below the simulated levels and store in the simulated array

/*          
            delta = tau[ich][ *isl-1] / tran[ *isl-1]  ;
            for(int il=0; il < *isl-1; ++il ) {
                tau[ich][il] = tran[il] * delta;
            }
            delta = tau[ich][ *iel-1] / tran[ *iel-1]  ;
            for(int il= *iel ; il < *ncel; ++il ) {
                tau[ich][il] = tran[il] * delta;
            }

*/

			delta = tau[ich][ *isl-1] - tran[ *isl-1]  ;
			for(int il=0; il < *isl-1; ++il ) {
				tau[ich][il] = tran[il] + delta;
			}
			delta = tau[ich][ *iel-1] - tran[ *iel-1]  ;
			for(int il= *iel ; il < *ncel; ++il ) {
				tau[ich][il] = tran[il] + delta;
			}

// save these transmissions
			old_tran.push_back( std::vector<double>(&tau[ich][0], &tau[ich][0]+ *ncel) );
/*
        for(int i=0; i< *ncel; ++i) {
          std::cout << ConvertToStringPrec(1.0-old_tran[ich][i]) << "   " << zaDouble[i] << std::endl;
        }
        exit(23);
*/
		}
        return;
	}
//    if(old_tran.empty() ) return;



	if( *bottom_up ) {    // bottom up  retrieval
		for( int ich = 0; ich< *nch; ++ich ) {
			delta = tau[ich][ *isl-1] - old_tran[ich][ *isl-1]  ;
			for(int il=0; il < *isl-1; ++il ) {
				tau[ich][il] = old_tran[ich][il] + delta;
			}

			if( *iel == *layer_stop) {   // need to scale below
				delta = tau[ich][ *iel-1] - old_tran[ich][ *iel-1]  ;
				for(int il= *iel ; il < *ncel; ++il ) {
					tau[ich][il] = old_tran[ich][il] + delta;
				}
			}

			old_tran[ich] = std::vector<double>(&tau[ich][0], &tau[ich][0]+ *ncel);
		}
	} else {     // top down retrieval

		for( int ich = 0; ich< *nch; ++ich ) {

			delta = tau[ich][ *iel-1] - old_tran[ich][ *iel-1]  ;
			for(int il= *iel ; il < *ncel; ++il ) {
				tau[ich][il] = old_tran[ich][il] + delta;
			}

			if( *isl == *layer_start) {   // need to scale above
				delta = tau[ich][ *isl-1] - old_tran[ich][ *isl-1]  ;
				for(int il=0; il < *isl-1; ++il ) {
					tau[ich][il] = old_tran[ich][il]  + delta;
				}
			}
			
			old_tran[ich] = std::vector<double>(&tau[ich][0], &tau[ich][0]+ *ncel);
		}

	}



/*
if( *iel - *isl > 10) {
	for(int il = 0; il < *ncel;  ++il ) {
		std::cout << *isl << "  " << *iel << " " << old_tran[0][il] << "  " << tau[0][il] << "  " << za[il] << std::endl;
	}
	exit(24);
}
*/

}