#include "F77Radtran.h"
#include "fortranfunctions.h"

#include <vector>
#include <numeric>
#include <iostream>
#include <cmath>
#include <cassert>
#include <algorithm>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>

struct function_params {
	int isl, iel;
	std::vector<float> za;
	std::vector<float> measured;
	std::vector<float> simulated;
};


extern "C" double bend_angle_rms(double delta, void *params)
{
//static int c__1 = 1;
	function_params* p = (function_params* ) params;

//int ncel = (int)p->za.size();
//std::cout << p->isl << "  " << p->iel << std::endl;
//    std::vector<double> diffs;
double sum = 0.0;
	for(int il=p->isl; il<= p->iel; ++il) 
		sum += std::pow( p->simulated[il]-(p->measured[il]+delta),2.0 );

    return  std::sqrt(sum / (p->iel - p->isl+1)  );  
/*
	double Mean = std::accumulate(diffs.begin(), diffs.end(), 0.0)/diffs.size() ; 

	double var=0.0;
    for(int i = 0; i<(int)diffs.size(); ++i) 
		var += std::pow( diffs[i] , 2.0);

std::cout << "delta " << delta*3600*180/M_PI << " mean "  << Mean*3600*180/M_PI << " std dev " 
          << std::sqrt(var / diffs.size() )*3600*180/M_PI  << std::endl;
	return std::sqrt(var / diffs.size() );
*/
}


extern "C" void create_merged_refrac(float* za,  float* sim_ref,
	float* meas_ref, int* ncel, float* zr, float* pr, float* tr, 
	float **tgr, int* nr, int* mr, float* re, float* merged_ref ,
        float* p_init, float* t_init, int* i_abel, int* i_abel2, 
	float* window_top, float* window_bottom)
/* sim ref    simulated refraction angle  (arcmin)
   meas_ref   measured bending angle   at altitudes za (arcmin)
   ncel       number of za points
  zreg   not used.
*/     
{

	static int c__7 = 7;
        static int c__1 = 1;
	static float m_off  = 5.0;

//std::cout << *nr << "  " << *ncel << "  " << za[0] << "  " << zr[0] << std::endl;
//exit(23);

// convert to arcseconds store in s and m
	std::vector<float> s(sim_ref, sim_ref+ *ncel), m( *ncel);
    std::transform(meas_ref, meas_ref+ *ncel , m.begin(), std::bind2nd( std::plus<float>(),m_off*M_PI/6.48e5 )  );


	std::vector<float>::reverse_iterator rit = std::find_if( s.rbegin(), s.rend(), std::bind2nd(std::less<float>(), (*window_bottom)*M_PI/6.48e5)) -1;
	int iel = std::distance(s.begin(), rit.base() );  // index of element at is >= window_bottom

	rit = std::find_if( s.rbegin(), s.rend(), std::bind2nd(std::less<float>(), (*window_top)*M_PI/6.48e5)) ;
	int isl = std::distance(s.begin(), rit.base() );  // index of element at is >=  window_top

//std::cout << " isl iel " << isl << " " << iel << std::endl;
//std::cout <<  *window_bottom << "  " << *window_top << std::endl; exit(23);

	std::vector<float> m_smooth(*ncel), s_smooth(*ncel), scale(*ncel);
	::bend_fit__( ncel, za, &m[0], &c__7, &c__7, za, &za[*ncel-1], &m_smooth[0], &scale[0] );
	::bend_fit__( ncel, za, &s[0], &c__7, &c__7, za, &za[*ncel-1], &s_smooth[0], &scale[0] );
/*
for(int j=0; j< *ncel; ++j)
    std::cout << s[j]* 180.*3600./M_PI << "  " << za[j] << std::endl;
std::cout << "&" << std::endl;

 for(int j=0; j< *ncel; ++j)
    std::cout  << m[j]* 180.*3600./M_PI << "  " << za[j] <<  std::endl;
std::cout << "&" << std::endl;
exit(23);
*/

/*
std::cout << " isl iel " << isl << " " << iel << std::endl;
std::cout <<  *window_bottom << "  " << *window_top << std::endl; exit(23);
*/

// gsl brent minimizer
	int status, iter=0;
	const static int max_iter=50;
        

	double mm, a= -1e-4 , b= 1e-4, guess = 1.e-5;
        const gsl_min_fminimizer_type *T;
        gsl_min_fminimizer *solver;
	gsl_function F;

	function_params params = {isl,iel, std::vector<float>(za, za+ *ncel), m_smooth,s_smooth};
	
	F.params = &params;
	F.function = &bend_angle_rms;

        T = gsl_min_fminimizer_brent;
        solver = gsl_min_fminimizer_alloc (T);
	assert(solver);
        status = gsl_min_fminimizer_set (solver, &F, guess, a, b );
	assert(status != GSL_EINVAL);

       do
         {
           iter++;
           status = gsl_min_fminimizer_iterate (solver);
     
           mm = gsl_min_fminimizer_x_minimum (solver);
           a = gsl_min_fminimizer_x_lower (solver);
           b = gsl_min_fminimizer_x_upper (solver);
     
           status 
             = gsl_min_test_interval (a, b, 1.e-7, 0.0);

         }
       while (status == GSL_CONTINUE && iter < max_iter);
       gsl_min_fminimizer_free (solver);

//if(status == GSL_SUCCESS) {
  assert(status == GSL_SUCCESS );
      std::cout << "merge window " << za[iel] << " to " << za[isl] << std::endl;
//  std::cout << "converged " << mm * 6.48e5/M_PI   << " arcseconds" << std::endl;
       params.measured = m;
       params.simulated = s;
// = {isl,iel, std::vector<float>(za, za+ *ncel), m,s};
      std::cout << "RMS is " << bend_angle_rms(mm, (void *)&params)*6.48e5/M_PI   << std::endl;
		mm += (m_off * M_PI/6.48e5) ;
		std::cout << "converged " << mm * 6.48e5/M_PI   << " arcseconds" << std::endl;
  assert( std::fabs(mm * 6.48e5/M_PI)  <= 0.5   ) ;
//      if (std::fabs(mm * 6.48e5/M_PI  ) > 0.5) 
//      {	std::cout << "Measured and Simulated refraction angles differ more than 0.5 arcseconds" << std::endl;  exit(24); }


//    } else {
//	   std::cout << "not converged" << std::endl;
//       exit(24);
//    }


/*   scale simulated
	float ratio = 1.0 - mm / s[iel];
	std::transform(s.begin(), s.begin()+iel+1 , merged_ref, std::bind2nd( std::multiplies<float>(),ratio )  );
	std::copy(m.begin()+iel+1, m.end(), merged_ref+iel+1);
*/

//  offset measured....

   std::transform(meas_ref, meas_ref+ *ncel, merged_ref, 
		std::bind2nd( std::plus<float>(),mm ) );


	std::copy(s.begin(), s.begin()+isl, merged_ref);


	float t = iel -isl ;
	for(int il = isl; il <= iel; ++il ) {
		float b = (il - isl)/ t;
		float a = 1.0 - b;
         merged_ref[il] = b * merged_ref[il] + a * s[il];

//		merged_ref[il] = a * merged_ref[il] + b * m[il];     //   scale simulated
	}

float *ll = std::find_if(&merged_ref[0], 
    &merged_ref[0]+ *ncel, std::bind2nd( std::greater_equal<float>(),10.0*M_PI/(180.*3600.) )  );
assert( ll != &merged_ref[0]+ *ncel );

ll = std::adjacent_find(ll, &merged_ref[0]+ *ncel,
	std::greater_equal<float>() );

if( ll != &merged_ref[0]+ *ncel) { 
	*ncel = (int)std::distance(&merged_ref[0] , ll )+1 ;
}

*i_abel = isl;
*i_abel2 = iel;


// to make the code write the entire profile out comment out the next line..


//std::cerr <<  za[0] << " " << za[*ncel-1] <<  " " << za [ (iel+isl)/2] << std::endl;
//exit(23);
//*zreg = za[iel];
//top_dz[0] = za[isl] - za[iel];
//*ztop = 
//if( *ztop <= 0) *ztop= za[ (isl+iel)/2] ;
//exit(23);
//*nsr = 1; *ner= *ncel;

/*
	float *pz = std::find_if(za,za+ *ncel, std::bind2nd(std::less_equal<float>(), zstart[0][0])) ;
	*nsr = std::distance(za, pz)+1;
	pz = std::find_if(za,za+ *ncel, std::bind2nd(std::less<float>(), zend[0][0]))-1 ;
	*ner = std::distance(za, pz)+1;
	if( *ner > *ncel) *ner = *ncel;
*/
//  find the pressure and temperature to use with this new merged bending angle
// We'll let path do the work


	static int debug = 0;   		// 0 is a false (this is the debug flag)
	static int irf =  0;    // don't need to do refraction here
	static int iplat = 0;  		// observer in space
//static int ic1 = -1;
//	static int ic2 = 1;
    float tearth,wave=1000., z2, z1 = za[ 0]; //+mm ;
//std::cout << "this is z1 " << z1 << std::endl;
	int mcel=  c__1  ;
    	int it2 = mcel * mcel;
    	int it1 = it2 + mcel;
    	int mpout = mcel * 11 + 3 + it1 * 12 + it2 * 3;
	std::vector<float> pout(mpout),gr(2*mcel);
	
	::kp_inita__(&debug, &mcel, (int*)&pout[0]  )  ;
	
//std::cout << " right here " << std::endl;
	::kp_mpath__( mr,  &mcel,  &mpout,  &irf,  &c__1, &c__1, &c__1, 
		re, &c__1, nr,
		&wave, &iplat, &z1, zr, pr, tr, &tgr[0][0],
 		&z2, p_init, t_init, &gr[0], &tearth, &pout[0] );

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

//exit(23) ; (*p_init) *= 0.2;
//std::cout <<  " mmmm " << std::endl; exit(23);
//std::cout << "tr zr " << std::endl;
//for(int l=0; l< *nr; ++l)
//  std::cout << tr[l] << "  " << zr[l] << std::endl;
//std::cout << "&" << std::endl;


/*

	std::cout << za[iel] << "  " << m[iel] << " " << m_smooth[iel] << " " << m_scale[iel] <<
        " " << s[iel] << " " << s_smooth[iel] << "  " << s_scale[iel] <<  std::endl;

exit(23);
*/
/*

	float *pz = std::find_if(za, za+ *ncel, std::bind2nd(std::less_equal<float>(), 73.));// 55.) );
	int iz = std::distance(za, pz);


	float ratio = meas_ref[iz]/ sim_ref[iz];
	std::transform(sim_ref, sim_ref+iz+1, merged_ref, std::bind2nd( std::multiplies<float>(),ratio )  );

	std::copy(meas_ref+iz+1, meas_ref+ *ncel, merged_ref+iz+1);
*/
/*
for(int j=0; j< *nr; ++j) 
	std::cout << tr[j] << "  " << zr[j] << std::endl;
std::cout << "&" << std::endl;
*/

/*
for(int j=0; j< *ncel; ++j) 
	std::cout << sim_ref[j]* 180.*3600./M_PI << "  " << za[j] << std::endl;
std::cout << "&" << std::endl;

 for(int j=0; j< *ncel; ++j) 
	std::cout  << (mm+meas_ref[j])* 180.*3600./M_PI << "  " << za[j] << std::endl;
std::cout << "&" << std::endl;
for(int j=0; j< *ncel; ++j) 
	std::cout  << merged_ref[j]* 180.*3600./M_PI  << "  " << za[j] << std::endl;
std::cout << "&" << std::endl;
exit(23);
*/

	return;
}