//  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'.

// -----------------------------------------------------------
//> US76 atmospheric model
//>
//> Temperature, pressure or density as function of altitude 
//> (altitude must be between 0 and 1000 km). 
//
// -----------------------------------------------------------

// Defaut values
hmin = 0.e3;
hmax = 1000.e3;
iplot = 3; 

desc_param = list(..
   CL_defParam("Altitude - min", hmin, units=['m', 'km'], id='$hmin', valid='$hmin>=0' ),..
   CL_defParam("Altitude - max", hmax, units=['m', 'km'], id='$hmax', valid='$hmax>$hmin & $hmax<=1000' ),..
   CL_defParam("Plot: 1=Temperature, 2=Pressure, 3=Density", iplot, accv=1:3)..
);

[hmin, hmax, iplot] = CL_inputParam(desc_param);


// -----------------------------------------------------------
// computation
// -----------------------------------------------------------

nbpts = 200;
alt = linspace(hmin, hmax, nbpts);

[temp, pres, dens] = CL_mod_atmUS76(alt);
res = [temp; pres; dens];


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------
// titles
texts = ["Temperature" , "Pressure" , "Density" ];
units = ["(K)" , "(Pa)" , "(kg/m^3)" ];


f = scf();
f.visible = "on";
f.immediate_drawing = "off";

y = res(iplot,:); 

plot(alt/1000, y, "thickness", 2);

// automatic log scale 
logscale = %f; 
if (min(y) > 0)
  if (max(y) / min(y) > 100); logscale = %t; end
end

a = gca(); 

if (logscale); a.log_flags = "nln"; end  

a.title.text = "US76 atmospheric model" + " - " + texts(iplot);
a.x_label.text = "Altitude (km)";
a.y_label.text = texts(iplot) + " " + units(iplot);
CL_g_stdaxes(a);

f.immediate_drawing = "on";
f.visible = "on";