#include "fortranfunctions.h"

#include <vector>
#include <algorithm>
#include <iostream>


void interp_2pt(float* xa, float* ya, float* y2, float* x2, int n1, int n2, float miss) 
{

/*
;------------------------------------------------------------------
;
; Do a 2 point linear interpolation of x,  from the y grid to the 
; y2 point (or grid) 
;
; input:
;   x......x axis variables 
;   y......y axis variables 
;   y2.....interpolate to this value or grid
;   miss...missing data value
;
; output:
;   result...the data value from x at y2
;
; source: Mark Hervig 
;------------------------------------------------------------------

  if (n_elements(x) ne n_elements(y)) then begin
    print,'x and y arrays are not the same size in interp_2pt.pro!!!'
    stop
  endif

  n  = n_elements(y)
  n2 = n_elements(y2)
  x2 = fltarr(n2) + miss

  ok = where(x ne miss,nok)

  std::fill(&x2[0], &x2[n2], miss);
   
  if (nok gt 1) then begin

    xa = x(ok)
    ya = y(ok)
*/
     float term;
    for(int i = 0; i <= n2-1; i++) { // do begin

      for (int j = 0; j<= n1-2; ++j) { //  do begin

        if (ya[j] >= y2[i] && ya[j+1] <= y2[i]) {// then begin
          term  = (y2[i]-ya[j]) / (ya[j+1]-ya[j]);
          x2[i] = term * (xa[j+1]-xa[j]) + xa[j];
          break; //goto,jump1
        }//endif
        if (ya[j] <= y2[i] && ya[j+1] >= y2[i]){ //then begin
          term  = (y2[i]-ya[j]) / (ya[j+1]-ya[j]) ;
          x2[i] = term * (xa[j+1]-xa[j]) + xa[j] ;
          break; //goto,jump1
        }//endif

       }//endfor

      
    } // endfor
    
}
/*
  endif

  return,x2
  end
*/



float ll_poly(const int n, const float* a, float x0){
   int i;
   a--;
   float p = a[n];
   for(i = n - 1; i >= 0; i--){
     p = a[i] + x0*p;
   }
   return p;
}


extern "C" void sofie_ice_t ( int* np, float* za, float* ext9, float* ext10, float* ext8,  float* ztop,
   float* zbot  )
{
   
   static std::vector<float> arx;
   
   static  int nar = 30, c__1=1, c__4=4 ;
   static  float Tm[] = {120.,130.,140.,150.}  ; //model temperatures (K)


   static  float fc_ar_r9_8[4][3] = {
    {  3.34158e+00, -1.50551e-01,  8.54202e-03}, // 120K
    {  3.03287e+00, -1.02951e-01,  5.83847e-03}, //130K
    {  2.72417e+00, -5.53503e-02,  3.13490e-03}, // 140K
    {  2.51530e+00, -3.02189e-02,  1.71940e-03}   } ; //150K
  
   static float fc_ar_r10_8[4][3] = {
   {  1.37762e+00,  2.07424e-01, -7.42188e-03} , // 120K
   {  1.13493e+00,  1.75585e-01, -5.94918e-03} ,// 130K
   {  8.92238e-01,  1.43745e-01, -4.47648e-03} ,// 140K
   {  8.15695e-01,  1.34612e-01, -4.07685e-03}  };  //150K
 
   static float fc_ar_r9_10[4][3] = {
   {  2.18975e+00, -3.13427e-01,  2.17300e-02} ,   // 120K
   {  2.55666e+00, -3.43872e-01,  2.33435e-02} , // 130K
   {  2.92357e+00, -3.74317e-01,  2.49571e-02} , //140K
   {  2.95423e+00, -3.66487e-01,  2.42033e-02}  }; // 150K
    

   if(arx.empty() ) {
       float dar = 5./(nar-1);
       for (int j=0; j< nar; ++j) 
	   arx.push_back( j*dar + 1.);
   }

for(int j =0; j< *np; ++j) {
     if (! (za[j] <= *ztop && ext9[j] > 5.e-7 && ext10[j] > 5.e-7  && ext8[j] > 5.e-7 )) continue;
     if( za[j] < *zbot ) break;

     float r9_8  = ext9[j] / ext8[j];
     float r10_8 = ext10[j] / ext8[j];
     float r9_10 = ext9[j] / ext10[j];
   
     std::vector<float> t9_8(nar,0.);
     std::vector<float> t10_8 = t9_8;
     std::vector<float> t9_10 = t9_8;
       
     for(int ir=0; ir< nar; ++ir) {
       float AR_in = arx[ir];
       std::vector<float> r9_8m, r9_10m, r10_8m;
     
       for(int ik = 0; ik <4; ++ik) {
	       r9_8m.push_back(  ll_poly(3,  &fc_ar_r9_8[ik][0],  AR_in )  );
	       r9_10m.push_back( ll_poly(3,  &fc_ar_r9_10[ik][0], AR_in )  );
	       r10_8m.push_back( ll_poly(3,  &fc_ar_r10_8[ik][0], AR_in )  );
	    }
     
        if (r9_8  >= *std::min_element(r9_8m.begin() , r9_8m.end() ) 
		  && r9_8  <= *std::max_element(r9_8m.begin(), r9_8m.end() ) )
	    	::liner_(  &r9_8m[0],&Tm[0], &r9_8, &t9_8[ir], &c__4, &c__1   );
	    	        
        if (r10_8 >= *std::min_element(r10_8m.begin() , r10_8m.end()  ) 
          && r10_8 <= *std::max_element(r10_8m.begin() , r10_8m.end() ) )
	        ::liner_(  &r10_8m[0],&Tm[0], &r10_8, &t10_8[ir], &c__4, &c__1   );
		         
        if (r9_10 >= *std::min_element(r9_10m.begin() , r9_10m.end()  ) 
          && r9_10 <= *std::max_element(r9_10m.begin() , r9_10m.end() ) ) 
	        ::liner_( &r9_10m[0], &Tm[0], &r9_10, &t9_10[ir], &c__4, &c__1   );
     }

// now find a common solution  Axial ratio vs Temp

	 std::vector<float> dt, ar, t9_8f;
	 for(int ir=0; ir< nar; ++ir) {
	    if (t9_8[ir] > 100. && t10_8[ir] > 100. ) {
	    	dt.push_back( t9_8[ir] - t10_8[ir] );
		    ar.push_back(arx[ir] );
		    t9_8f.push_back(t9_8[ir]);
	    }
     }
	
//interp_2pt(float* xa, float* ya, float* y2, float* x2, int n1, int n2, float miss) 	
     float zero = 0.0, ar_ret, t_ret;
	 int nk = (int)dt.size();

//std::cout << za[j] << "  " << nk << std::endl;

	if(nk > 1) {
	 interp_2pt(  &ar[0], &dt[0], &zero, &ar_ret, nk, 1, -999. );
	 interp_2pt(  &t9_8f[0], &ar[0], &ar_ret, &t_ret, nk, 1, 0. );
	}
	
//;- Determine the band 7/8 & 7/9 extinction ratios for t_ret

    float r7_8m[] = {0.0195491,   0.0275128,    0.0354764,    0.0623039} ; //this is vs. Tm, for AR=2 (not a f(AR))
    float r7_9m[] = {0.0057458,   0.0096184,    0.0134910,    0.0252974} ; // "
   
    float r7_8 = 0. ;
    float r7_9 = 0.  ;
   
    if (t_ret >= 120. &&  t_ret <= 150.)  { // then begin    
      ::liner_(  &Tm[0],&r7_8m[0],&t_ret, &r7_8,  &c__4, &c__1   ); 
      ::liner_(  &Tm[0],&r7_9m[0],&t_ret, &r7_9,  &c__4, &c__1   ); 
	
      std::cout << za[j] << " " << t_ret << "  " << ar_ret << "  " << r7_8 << "  " << r7_9 
<< "   " << ext9[j] << "  " << ext10[j] << "  " << ext8[j]   << std::endl;
    }

} // end loop over layers

exit(23);


}