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

// -----------------------------------------------------------
//> Various visibility parameters.  
//> See CL_gm_visiParams for details. 
//
// Author: A. Lamy
// -----------------------------------------------------------

// Defaut values
tab_alt = [500, 700, 1000, 1400, 2000, 5000, 10000, 20000, 36000] * 1.e3;
angtype = 2; 
angmin = 0 * %CL_deg2rad;
angmax = 90 * %CL_deg2rad;
show = 3;
elevmark = 5 * %CL_deg2rad; 

desc_param = list(..
   CL_defParam("Levels of altitude", tab_alt, units=['m', 'km'], dim=[1,15], valid='$x>0' ),..
   CL_defParam("Input angle: 1=elevation, 2=center angle, 3=sat angle, 4=incidence", angtype, id="$angtype", accv=1:4),..
   CL_defParam("Input angle value - min", angmin, units=['rad', 'deg'], id='$angmin', valid='$angmin >= 0 & $angmin <= 90' ),..
   CL_defParam("Input angle value - max", angmax, units=['rad', 'deg'], id='$angmax', valid='$angmax <= 90 & $angmax > $angmin' ),..
   CL_defParam("Plot parameter: 1=elevation, 2=center angle, 3=sat angle, 4=incidence, 5=distance", show, accv=1:5, valid="$x <> $angtype"),..
   CL_defParam("Elevation for markers (-1 => no marker)", elevmark, units=['rad', 'deg'], valid='($x >= 0 & $x <= 90) | $x == -1')..
);

[tab_alt,angtype,angmin,angmax,show,elevmark] = CL_inputParam(desc_param);

// -----------------------------------------------------------
// computation
// -----------------------------------------------------------
NAME = ["elev", "cen", "sat", "incid", "dist"]; 
ID = 1:5; 

nbpts = 1000; 
ang = linspace(angmin, angmax, nbpts); 

result = list(); 

for k = 1 : length(tab_alt)
  alt = tab_alt(k); 
 
  // computation for elevation [0, 0.5, 1] for a better accuracy 
  elev0 = [0 : 0.1 : 5] * %CL_deg2rad; 
  ang0 = CL_gm_visiParams(%CL_eqRad + alt, %CL_eqRad, "elev", elev0, [NAME(angtype)]);
  
  angres = ang; 
  I = find(ang0 >= angmin & ang0 <= angmax); 
  angres = CL_sortMat([ang, ang0(I)]); 
    
  [res, elev] = CL_gm_visiParams(%CL_eqRad + alt, %CL_eqRad, NAME(angtype), angres, [NAME(show), "elev"]);
  I = find(elev < 0 | isnan(elev)); 
  res(I) = %nan; 
  angres(I) = %nan; 
  
  // computation for elevation "elevmark"
  [angm, resm] = CL_gm_visiParams(%CL_eqRad + alt, %CL_eqRad, "elev", elevmark, [NAME(angtype), NAME(show)]);
  if (angm < angmin | angm > angmax | elevmark < 0)
    angm = %nan; 
    resm = %nan; 
  end
    
  result(k) = struct("ang", angres, "res", res, "angm", angm, "resm", resm); 
end


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------
UNIT = [%CL_rad2deg, %CL_rad2deg, %CL_rad2deg , %CL_rad2deg, 1/1000];
UNIT_STR = ["deg", "deg", "deg", "deg", "km"]; 

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

f.color_map = jetcolormap(max(length(tab_alt), 5));

for k = 1 : size(result);
  plot(result(k).ang * UNIT(angtype), result(k).res * UNIT(show), "thickness", 2);
  h = CL_g_select(gce(), "Polyline"); 
  h.foreground = k; 
end

for k = 1 : size(result);
  plot(result(k).angm * UNIT(angtype), result(k).resm * UNIT(show), "o");
  h = CL_g_select(gce(), "Polyline"); 
  h.mark_foreground = color("black"); 
  h.mark_background = k; 
end


// general setting
titles = ["Elevation", "View angle from Earth centre", "View angle from satellite", "Incidence", "Distance"];

a = gca(); 

CL_g_stdaxes(a, colg=color("grey60"));

resmin = %inf; 
resmax = -%inf; 
angmin = %inf; 
angmax = -%inf; 

for k = 1 : size(result);
  resmin = min(resmin, min(result(k).res)); 
  resmax = max(resmax, max(result(k).res)); 
  angmin = min(angmin, min(result(k).ang));
  angmax = max(angmax, max(result(k).ang)); 
end

if (~isinf(angmin) & ~isinf(angmax))
  a.data_bounds = [angmin*UNIT(angtype),resmin*UNIT(show);angmax*UNIT(angtype),resmax*UNIT(show)];
  a.tight_limits="on";
end

a.title.text = titles(show);
a.x_label.text = titles(angtype) + " (" + UNIT_STR(angtype) + ")";
a.y_label.text = titles(show) + " (" + UNIT_STR(show) + ")";

CL_g_legend(a, string(tab_alt/1000), header="Altitude (km)"); 

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