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


h = 700.e3; // altitude
elevmin = 0 * %CL_deg2rad;
tab_elevmax = [10:10:90] * %CL_deg2rad;
x_axis = 1; 

desc_param = list(..
   CL_defParam("Altitude", h, units=['m', 'km'], id='$h', valid='$h>=0' ),..
   CL_defParam("Minimum elevation for visibility", elevmin, units=['rad', 'deg'], id='$elevmin', valid='$elevmin>=0 & $elevmin<=90' ),..
   CL_defParam("Levels of maximum elevation", tab_elevmax, units=['rad', 'deg'], dim=[1, 20], id='$elevmax', valid='$elevmax>=$elevmin & $elevmax<=90'),..
   CL_defParam("x-axis (1=time, 2=centre angle", x_axis, accv=[1,2])..
);

[h,elevmin,tab_elevmax,x_axis] = CL_inputParam(desc_param);

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

sma = h + %CL_eqRad; 
T = CL_kp_params('per', sma); 
nbpts = 200; // nombre de points en x

for k = 1:length(tab_elevmax);

  elevmaxpass = tab_elevmax(k); 
  tab_elev = linspace(elevmin,elevmaxpass,nbpts);

  angc1 = ang_visi(sma, tab_elev, elevmaxpass);
  t1 = -(angc1/(2*%pi)) .* T / 60;    // -> minutes

  t = [t1, -t1];
  angc = [-angc1, angc1];
  elev = [tab_elev, tab_elev]; 

  [t, I] = gsort(t, 'c', 'i'); 
  elev = elev(I);  
  angc = angc(I);  

  if (x_axis == 1)
     plot2d(t, elev*%CL_rad2deg, style=k);
  else
     plot2d(angc*%CL_rad2deg, elev*%CL_rad2deg, style=k);
  end
end


// 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 pass - non rotating Earth (alt = " + string(h/1000)+ " km)";
if (x_axis == 1)
  a.x_label.text = "Time (mn)";
else
  a.x_label.text = "Centre angle (deg)";
end
a.y_label.text = "Elevation (deg)";

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

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