#include "Level1Data.h"
#include "fortranfunctions.h"

#include <cassert>
#include <iostream>
#include <cmath>


extern "C" void get_meas_refraction_angle( float*  za, float* refrac_angle,
	int* ncel, float **zstart, float **zend, int* nchn,	int* ngdo,
	int* nsr, int* ner)
{

std::cout << " inside refrac " << std::endl;

//std::cout <<  za[0] << "  " << za[*ncel-1] << std::endl; exit(23);
	
	std::vector<double> altgrid =  ::gLevel1Data.GetImpactAlts() ;
	std::vector<double> refangle = ::gLevel1Data.SunSensorRefraction();

	assert(!refangle.empty() && altgrid.size() == refangle.size() );

	std::vector<float> x,y;

	float z1 = -1000., z2=1000.;
	for(int i = 0; i < *ngdo; ++i) {
		for(int j = 0; j < *nchn; ++j) {
			z1 = std::max(z1, zstart[j][i]);
			z2 = std::min(z2, zend[j][i]);
		}
	}

//	std::cout << *ngdo << "   " << *nchn << "  " <<  z1 << "  "  << z2 << std::endl ;
//	exit(23);

	for(unsigned i = 0; i< altgrid.size(); ++i) {
	  //	  std::cout << altgrid[i] << "  " << refangle[i] << std::endl;
		if(//altgrid[i] >= z2 && altgrid[i] <= z1 && 
                refangle[i] > -1. ) {
			x.push_back(refangle[i] );
			y.push_back(altgrid[i]);
		}
	}
        assert(!x.empty());
//	altgrid=y;
//	refangle=x;
	int num1 = x.size();
 /*     for(int i = 0; i < num1; ++i)
        std::cout << y[i] << "  " << x[i] << std::endl;
   exit(23);
*/
//std::vector<double> zaDouble(za, za+ *ncel);


	float *pz = std::find_if(za, za+ *ncel,
		std::bind2nd(std::less_equal<float>(),  std::min(z1, y.front()) ) ) ;

	float *pb = std::find_if(za, za+ *ncel,
		std::bind2nd(std::less<float>(),  std::max(z2, y.back()) ) )-1 ;

//   std::cout << *pb  << std::endl ; exit(23);
	std::copy(pz, pb+1, za);

//        float *pz = std::find_if(za, za+ *ncel,
//		std::bind2nd(std::less<float>(),  y.back()  ) ) ;

        *ncel = (int)std::distance(pz, pb)+1;
// std::cout << zaDouble[ *ncel-1] << "   " << altgrid.back() << std::endl; 

 	::liner_(&y[0], &x[0], za, refrac_angle, &num1, ncel );
    *nsr= 1;
	*ner = *ncel;

 std::cout << "nsr ner " << *nsr << "  " << *ner << std::endl;
//      for(int i = 0; i < *ncel; ++i) 
//        std::cout << za[i] << "  " << refrac_angle[i] << std::endl;
//      exit(23);

}