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

// -----------------------------------------------------------
//> Visibility duration:
//> Duration of a visibility pass from a user on the ground. 
//> NB: approximate (geometrical) calculation, valid only for low orbits. 
//> Hypothesis: spherical and non rotating Earth
//
// Author: A. Lamy
// -----------------------------------------------------------


hmin = 400.e3;
hmax = 1400.e3;
elevmin = 5 * %CL_deg2rad;
tab_elevmax = [7.5, 10, 15, 20, 30, 50, 90] * %CL_deg2rad;

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' ),..
   CL_defParam("Minimum elevation for visibility", elevmin, units=['rad', 'deg'], id='$elevmin', valid='$elevmin>=0' ),..
   CL_defParam("Levels of maximum elevation", tab_elevmax, units=['rad', 'deg'], dim=[1, 20], id='$elevmax', valid='$elevmax>=$elevmin & $elevmax <=90' )..
);

[hmin,hmax,elevmin,tab_elevmax] = CL_inputParam(desc_param);

nbpts = 60; // nombre de points en x
tab_h = linspace(hmin,hmax,nbpts);

// -----------------------------------------------------------
// half visibility angle
// elev_min : elevation at the beginning/end of each pass
// elev_max_pass : max elevation of the pass
// NB : could be optimized by pre-compiting a few data
// (not done to keep everything in one function)
// -----------------------------------------------------------

function [ang] = ang_visi(sma, elevmin, elevmaxpass)
  N = max(length(sma), length(elevmin), length(elevmaxpass)); 

  sma = sma .* ones(1,N);
  elevmin = elevmin .* ones(1,N);
  elevmaxpass = elevmaxpass .* ones(1,N);

  r0 = %CL_eqRad * ones(1,N); 

  // angle from Earth center at beginning of pass
  cen_beg = CL_gm_visiParams(sma, r0, 'elev', elevmin, 'cen'); 
  // angle from Earth center at mid-pass
  cen_mid = CL_gm_visiParams(sma, r0, 'elev', elevmaxpass, 'cen'); 

  ang = acos(cos(cen_beg) ./ cos(cen_mid)); 
endfunction

// -----------------------------------------------------------
// results / plot
// -----------------------------------------------------------

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

a=gca();

nb = length(tab_elevmax);
f.color_map = jetcolormap(min(max(nb,3),100)) * 0.95; // min 3 valeurs

tab_sma = tab_h + %CL_eqRad; 
T = CL_kp_params('per', tab_sma); 

duree_min = %inf; 
duree_max = -%inf; 

for k = 1:length(tab_elevmax);

  elevmax_pass = tab_elevmax(k);

  ang = ang_visi(tab_sma, elevmin, elevmax_pass);
  duree = 2*(ang/(2*%pi)) .* T / 60;    // -> minutes

  duree_min = min(duree_min, min(duree));
  duree_max = max(duree_max, max(duree));

  plot2d(tab_h/1000, duree , style=k);
end


// legends
a1 = CL_g_legend(a, string(tab_elevmax*%CL_rad2deg) + " deg", header="Max elevation :");


// general setting
CL_g_stdaxes(a)
a.data_bounds = [hmin/1000,duree_min; hmax/1000,duree_max];
a.tight_limits="on";
a.title.text = "Visibility duration - non rotating Earth (min elev = " + string(elevmin*%CL_rad2deg)+ " deg)";
a.x_label.text = "Altitude (km)";
a.y_label.text = "Duration (mn)";

// change properties
h=CL_g_select(a, "Polyline");
h.thickness=2;

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