#include "frost_point.h"

#include <algorithm>
#include <cassert>
#include <cmath>

void ice_equilib3 (float* p_in, float* t_in, float* q_h2o,
	int* opt, float* t_ice, float* p_ice, float* v_ice,  int* nt)
{
/*
;-------------------------------------------------------
;  purpose:
;
;     Routine computes the frost point temperature, saturation vapor 
;     pressure over ice at T_in, and volume of condensed ice.
;    
;
;  input:
;     t_in....temp, K
;     p_in....pressure, mb
;     q_h2o...h2o mixing ratio, ppmv
;     opt.....1= use Rapp 2005, 2= use Mauersberger 2003
;     nt   number of elements
;
;  output:
;     t_ice...frost point temperature, K
;     p_ice...saturation vapor pressure over ice at T_in, mb
;     v_ice...volume of condensed ice, um3/cm3
;
;  requires goff_gratch_ice.pro
;
;  source: Mark Hervig 
;-------------------------------------------------------------


;- setup some variables
*/

   const static float Mww   = 18.0                 ; //molec wt. h2o, g/mol
   const static float di    = 0.916                ; //density ice, g/cm3
   const static float R     = 8.314                ; //J/mol/K
/*
   nt = n_elements(t_in)
   np = n_elements(p_in)
   nq = n_elements(q_h2o)

   if (nt ne np or nt ne nq or nq ne np) then begin
     print,'temperature,  pressure,  and h2o must have the same dimensions in ice_equilib_sphere.pro'
     stop
   endif
*/
  
   assert( *opt == 1 ||  *opt == 2); 
//   if (opt ne 1 and opt ne 2) then stop,'opt not valid in ice_equilib3.pro'

	std::fill(t_ice, t_ice+ *nt,   -1.0e+24);   // fill with missing values
	std::fill (p_ice, p_ice+ *nt,  -1.0e+24);
	std::fill (v_ice, v_ice+ *nt,  -1.0e+24);

//   t_ice = fltarr(nt) -999.
//   p_ice = fltarr(nt) -999.
//   v_ice = fltarr(nt)
   
//;- do the calculations

	float p1, p2;
//   for i = 0,nt-1 do begin
	for(int i=0; i< *nt; ++i) {
		if(t_in[i] >= 273.15) continue ;
//    if (t_in(i) lt 273.15) then begin

     		float p_h2o = p_in[i] * q_h2o[i] *1.e-6  ; //h2o partial pressure, mb
     		float dpmin = 0.07 * p_h2o               ; // find t_ice to within 7% of vap pres

/*
;-   find the frost point temperature t_ice.  Use the fact that
;    a linear relationship exists between 1/T and alog10(p_h2o).
;    Do a linear interpolation between some hi and low temps
*/

     		float t1 = t_in[i] - 20. ; //173.15 
     		float t2 = t_in[i] + 20. ; //273.15 
          
     		if (*opt == 1) {  //then begin
       			p1 = std::log10( h2o_vapr_ice_rapp05(t1)  );
       			p2 = std::log10( h2o_vapr_ice_rapp05(t2)  );
     		} else {
       			p1 = std::log10( h2o_vapr_ice_mauersberger(t1)  );
       			p2 = std::log10( h2o_vapr_ice_mauersberger(t2)  );
    		}
     
     		float ta = 1./t1;
     		float tb = 1./t2;
     		float tc = 1./t_in[i];

     		float dsds     = (p2 - std::log10( p_h2o ) ) / (p2 - p1) ;
     		t_ice[i] = 1.0 / ( tb - dsds * (tb - ta) ) ;

//;-   get the saturation vapor pressure at t_in
          	
		p_ice[i] = (*opt == 1) ? h2o_vapr_ice_rapp05( t_in[i] ) : h2o_vapr_ice_mauersberger( t_in[i] ) ;
/*
     if (opt eq 1) then p_ice(i) = h2o_vapr_ice_rapp05( t_in(i) ) 
     if (opt eq 2) then p_ice(i) = h2o_vapr_ice_mauersberger( t_in(i) )         

;-   if saturated, find condendsed (equilibrium) ice volume
*/

     		float q_xs = q_h2o[i]*1.e-6 - p_ice[i] / p_in[i]   ; //excess h2o mix ratio
     		float n_xs = p_in[i] *100.0 * q_xs / (R*t_in[i] )   ; //excess mols h2o per m3 air

     		if (q_xs > 0.0) v_ice[i]  = 1.0e6 * (n_xs * Mww) / di ;
	}
  //  endif

  // endfor
   
//;- done

//   return
  // end
}