//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [app_elev] = CL_mod_atmRefract(true_elev, temp, pres, mod)
// Effect of atmospheric refraction on elevation from Earth surface
//
// Calling Sequence
// [app_elev] = CL_mod_atmRefract(true_elev [, temp, pres, mod])
//
// Description
// <itemizedlist><listitem>
// <p>Effect of refraction on elevation due to Earth atmosphere (for visible light) as seen from 
// Earth surface.</p>
// <p>The function computes the apparent elevation from the true elevation ("true" meaning:
// "that would be observed without atmosphere"). Because of refraction, a body actually below the horizon 
// can still be visible. </p>
// <p>Two models are available: Benett or Saemudsson. </p>
// <p>Effect of temperature and pressure can (roughly) be taken into account. </p>
// <p></p></listitem>
// <listitem>
// <p>Notes: </p>
// <p>- The refraction models come from Meeus (see Bibliography below). </p> 
// <p>- The refraction corrections are slightly adjusted so that there are exactly 0 
// for an elevation of 90 deg and unchanged for a true elevation of 0 deg. </p>
// <p>- Saemundson's formula directly gives the correction to be applied to the true elevation. But 
// Bennett's formula is to be applied to the apparent elevation. The formula is inversed
// so that the correction can be computed from the true elevation. 
// The effect of temperature and pressure is applied to both models, allthough Meeus
// mentions it for Saemundsson's model only. </p>
// <p>- The refraction corrections are computable for true elevations above -1 deg. 
// The correction is undefined (%nan) below this value. </p>
// <p> The 2 models are consistent to within less than 7 arcseconds between apparent elevations 
// 0 and 90 degrees. Benett model is supposed to be slightly more accurate (accuracy assumed better than 
// 4 arcseconds between 0 and 90 deg (apparent elevation). </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// true_elev: True elevation [rad]. (1xN or 1x1)
// temp: (optional) Atmospheric temperature [K]. Default is [] which means 283.15 K (10 °C). (1xN or 1x1)
// pres: (optional) Atmospheric pressure [Pa]. Default is [] which means 1.01e5 Pa (1010 mbar). (1xN or 1x1)
// mod: (string, optional) Model used: "ben" = Benett, "sae" = Saemudsson. Default is "ben". (1x1)
// app_elev: Apparent elevation [rad]. (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Astronomical algorithms, Jean Meeus, 2nd Edition, p 106  
// 
// Examples
// elev = (0 : 90) * (%pi/180); 
// r = CL_mod_atmRefract(elev) - elev; 
// scf(); 
// plot(elev * (180/%pi), r * (180/%pi), "b"); 
//

// Declarations: 

// Saemundsson model
// h: TRUE elevation (rad)
// r = refraction correction = apparent elevation minus true elevation (>= 0) (rad)
function [r] = refrac_saemundsson(h)
  c = %pi/180; 
  r = (1.02 ./ tan((h/c + 10.3 ./ (h/c + 5.11)) * c)) * (c / 60);
endfunction

// Saemundsson model + adjustment (0 for elev = 90 deg, unchanged for elev = 0 deg)
// h: TRUE elevation (rad)
// r = refraction correction = apparent elevation minus true elevation (>= 0) (rad)
function [r] = refrac_saemundssonC(h)
  h_all = [h, 0, %pi/2]; 
  r_all = refrac_saemundsson(h_all); 
  r = r_all(1:$-2); 
  r0 = r_all($-1); 
  r90 = r_all($); 
  r = (r - r90) * (r0 / (r0 - r90));
endfunction

// Benett model
// h: APPARENT elevation (rad)
// r = refraction correction = apparent elevation minus true elevation (>= 0) (rad)
function [r] = refrac_benett(h)
  c = %pi/180; 
  r = (1 ./ tan((h/c + 7.31 ./ (h/c + 4.4)) * c)); // arcmin
  // benett additional refinement
  dr = -0.06 * sin((14.7 * r + 13) * c); // arcmin
  r = (r + dr) * c / 60; 
endfunction

// Benett model inversed
// h: TRUE elevation (rad)
// r = refraction correction = apparent elevation minus true elevation (>= 0) (rad)
function [r] = refrac_benettInv(h)
  // The function solves: app_elev = true_elev + corr = true_elev + benett(app_elev)
  // => app_elev ? such that: F(app_elev) = app_elev - benett(app_elev) - true_elev = 0  
  function [z] = refrac_F(happ, I, args)
    // args: true elevation (1xN)
    z = happ - refrac_benett(happ) - args(I); 
  endfunction

  delta = 1 * %pi/180; // assumed maximum correction = 1 deg
  args = h; 
  // bounds (min/max)
  h_bounds1 = max(h-delta, -%pi/2); 
  h_bounds2 = min(h+delta,  %pi/2);
  happ = CL_fsolveb(refrac_F, h_bounds1, h_bounds2, args, ytol = 1.e-14, meth="ds"); 
  r = happ - h; 
endfunction

  // Benett inversed + adjustment (0 for elev = 90 deg, unchanged for elev = 0 deg) 
// h: TRUE elevation (rad)
// r = refraction correction = apparent elevation minus true elevation (>= 0) (rad)
function [r] = refrac_benettInvC(h)
  h_all = [h, 0, %pi/2]; 
  r_all = refrac_benettInv(h_all); 
  r = r_all(1:$-2); 
  r0 = r_all($-1); 
  r90 = r_all($); 
  r = (r - r90) * (r0 / (r0 - r90));
endfunction

// additional correction factor function of temperature and pressure
// T: temperature (K)
// P: pressure (Pa) 
function [f] = refrac_factTP(temp, pres)
  P = (pres / 1.e5) * 1000; // millibars
  T = temp - 273.15; // degrees Celcius
  f = (P / 1010) .* (283 ./ (273 + T)); 
endfunction
  
  
// Code:

// defaut values so that correcting factor is 1
if (~exists("temp","local")); temp = []; end
if (~exists("pres","local")); pres = []; end
if (~exists("mod","local")); mod = "ben"; end

// Defaut values chosen so that correction factor is 1
if (temp == [])
  temp = 283.15; // K
end

if (pres == [])
  pres = 1.01e5; // Pa
end

// Check model name
if (mod <> "ben" & mod <> "sae")
  CL__error("Invalid model"); 
end

// Check elevation value
if (find(true_elev < -%pi/2 - 10 * %eps | true_elev > %pi/2 + 10 * %eps) <> [])
  CL__error("Invalid value for argument elev"); 
end

// Check size / resize
[true_elev, temp, pres] = CL__checkInputs(true_elev, 1, temp, 1, pres, 1); 

if (true_elev == [])
  app_elev = []; 
  return; 
end

// minimum elevation value for model validity (arbitrary)  
elevb = -1 * %pi/180; 

I = find(true_elev < elevb)
true_elev(I) = %nan; 

if (mod == "ben")
  refrac = refrac_benettInvC(true_elev); 
else
  refrac = refrac_saemundssonC(true_elev); 
end

app_elev = true_elev + refrac .* refrac_factTP(temp, pres);

endfunction