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

// -----------------------------------------------------------
//> Local time at the satellite position 
//> (as function of latitude or argument of latitude)
//> The graph shows the difference between the satellite local time
//> and the mean local time of the ascending node (MLTAN). 
//> The min and max values correspond to the effect of the equation of time.    
//> NB: 
//> The x-axis in the case 'latitude - increasing' is reversed. 
//> (argument of latitude is increasing from left to right) 
//
// Author: A. Lamy
// -----------------------------------------------------------

// -----------------------------------------------------------
// Utility functions
// -----------------------------------------------------------
// inc: inclination (1x1)
// mltan: mean local time of ascending node (hours)
// x: latitude or arg of latitude 
// xopt: meaning of x: 1=lat-asc, 2=lat-desc, 3=pso-asc, 4=pso-desc
// t: time (1950.0) 
// NB: t or x is 1xN

function [res] = sat_info(inc, mltan, x, xopt, t)
  // names corresponding to xopt for CL_gm_orbSphTriangle
  xopt_vals = [1, 2, 3, 4]; 
  xopt_names = ["decl_a", "decl_d", "aol", "aol"]; 
  
  N = max(length(t), length(x)); 
  t = t .* ones(1,N); 
  x = x .* ones(1,N); 
    
  [pso, lat, alpha] = CL_gm_orbSphTriangle(inc, xopt_names(xopt), x, output=["aol", "decl", "rra"]);

  // RAAN (t is made 1xN) => output is 1xN
  gom = CL_op_locTime(t, 'mlh', mltan, 'ra'); 

  usat = CL_co_sph2car([gom+alpha; lat; ones(lat)]); 
  pos_sun = CL_eph_sun(t); 

  sza = CL_vectAngle(usat, pos_sun);  
  [tlh, mlh] = CL_op_locTime(t, 'ra', gom+alpha, ['tlh', 'mlh']); 

  res = struct(); 
  res.sza = sza;
  res.mlh = mlh; 
  res.tlh = tlh; 
  res.gom = gom; 
  res.pso = pso; 

endfunction


// -----------------------------------------------------------
// Parameters input
// -----------------------------------------------------------

inc = 98 * %CL_deg2rad; 
xopt = 1; 

desc_param = list(..
   CL_defParam("Orbit inclination", inc, units=['rad', 'deg'], id='$inc', valid='$inc > 0 & $inc < 180'),..
   CL_defParam("x-axis (1=lat-asc, 2=lat-desc, 3=pso-asc, 4=pso-desc)", xopt, accv=1:4)..
);

[inc, xopt] = CL_inputParam(desc_param);


// -----------------------------------------------------------
// computations
// -----------------------------------------------------------

if (xopt == 1 | xopt == 2)
  xmin = -min(inc, %pi-inc); 
  xmax = -xmin; 
elseif (xopt == 3) 
  xmin = -%pi/2; 
  xmax = %pi/2; 
else  
  xmin = %pi/2; 
  xmax = 3*%pi/2; 
end

tabx = linspace(xmin, xmax, 1000); 

hl = 0; // mean local time of asc node (arbitrary)
t0 = CL_dat_cal2cjd(2012,1,1); // arbitrary

res = sat_info(inc, hl, tabx, xopt, t0); 

if (xopt == 1 | xopt == 3)
  mlh = CL_rMod(res.mlh, -12, 12); 
else
  mlh = CL_rMod(res.mlh, 0, 24); 
end

// compute min/max true local time (equation of time)
t = t0:0.1:t0+365; 
dlt = CL_rMod(CL_op_locTime(t, 'mlh', 0, 'tlh'), -12, 12); 
tlhmin = CL_rMod(mlh+min(dlt), mlh-12, mlh+12); 
tlhmax = CL_rMod(mlh+max(dlt), mlh-12, mlh+12); 

function [val] = suppr_jump(val)
  I = find(abs(val(2:$) - val(1:$-1)) > 12); 
  val(I) = %nan; 
endfunction

mlh = suppr_jump(mlh); 
tlhmin = suppr_jump(tlhmin); 
tlhmax = suppr_jump(tlhmax); 


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------

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

plot(tabx*%CL_rad2deg, mlh, "k"); 
plot(tabx*%CL_rad2deg, tlhmin, "b"); 
plot(tabx*%CL_rad2deg, tlhmax, "r"); 


// plot adjustments
a=gca();

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

// titles
x_labels = ["Latitude (deg) - asc", "Latitude (deg) - desc", "pso (deg) - asc", "pso (deg) - desc"]; 
a.x_label.text = x_labels(xopt); 
a.y_label.text = "Local time minus MLTAN (h)"; 
a.title.text = sprintf("Satellite local time / MLTAN (h) - inc = %.1f deg", inc * %CL_rad2deg); 

// reverse x axis if latitude descending
// y_location has to be changed ! 
if (xopt == 2)
  a.axes_reverse(1) = "on"; 
end

a.data_bounds = [xmin*%CL_rad2deg, min(tlhmin)-0.1; xmax*%CL_rad2deg, max(tlhmax)+0.1]; 
a.tight_limits = "on"; 

CL_g_stdaxes(a);

CL_g_legend(a, str=["Mean", "True (min)", "True (max)"]); 

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