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

// -----------------------------------------------------------
//> Solar Zenith Angle: 
//> is defined (here) as the angle between the Satellite radius vector 
//> and the "Earth center-> Sun" vector (so that the altitude of the satellite 
//> has no impact). 
//>
//> NB: 
//> "MLTAN" means: mean local time of ascending node. 
//
// 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; 
hl = 10.5; // heure locale moyenne of asc node
tablat = (0:10:70) * %CL_deg2rad; 
asc_desc = 1; 

desc_param = list(..
   CL_defParam("Orbit inclination", inc, units=['rad', 'deg'], id='$inc', valid='$inc > 0 & $inc < 180'),..
   CL_defParam("Mean local time of ascending node", hl, units=['h'], valid='$x>=0 & $x<=24'),..
   CL_defParam("Latitudes", tablat, units=['rad', 'deg'], valid='abs($x) >= -90 & abs($x) <= 90', dim = [1, 20]),..
   CL_defParam("Ground tracks (1=ascending, 2=descending)", asc_desc, accv=[1,2])..
);

[inc, hl, tablat, asc_desc] = CL_inputParam(desc_param);


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

date_now = CL_dat_now("cal");
year = date_now(1);  // year number (current year)

t0 = CL_dat_cal2cjd(year,1,1); 
tf = CL_dat_cal2cjd(year+1,1,1); 
yearlen = tf-t0; 

s = linspace(0,1,50); // time in year (between 0 and 1)
t = t0 + s * yearlen ; 

// check tablat
I = find(abs(tablat) > min(%pi, %pi-inc)); 
tablat(I) = [];  

sza = zeros(length(tablat), length(t)); 

for k = 1:length(tablat);
  lat = tablat(k); 
  res = sat_info(inc, hl, lat, asc_desc, t); 
  sza(k,:) = res.sza; 
end


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

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

nbcol = max(length(tablat), 3); 
f.color_map = jetcolormap(nbcol); 

for k = 1:length(tablat);
  lat = tablat(k); 
  plot2d(s, sza(k,:) * %CL_rad2deg, style=k); 
end

// plot adjustments

a = gca(); 

// titles
a.x_label.text = "Day / Month (" + string(year) + ")"; 
a.y_label.text = "Solar zenith angle (deg)"; 

hl_cal = CL_dat_cjd2cal(hl/24); 
a.title.text = sprintf("Solar zenith angle - inc = %.1f deg, MLTAN = %dh%2d", inc*%CL_rad2deg, hl_cal(4), hl_cal(5)); 
CL_g_stdaxes(a);


// ticks
m = [1:12, 1]; 
d = ones(m);
y = [year * ones(1:12), year+1]; 
t_ticks = CL_dat_cal2cjd(y, m, d); 
labels = "1/" + string(m); 
x_ticks = (t_ticks - t0) / yearlen; 

a.x_ticks = tlist(['ticks', 'locations', 'labels'], x_ticks, labels);
 

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

CL_g_stdaxes(a);
a.sub_ticks(1) = 0; // no substicks for a-axis 
a.grid_position = "background"; 

str = ["asc", "desc"]; 
CL_g_legend(a, str = string(tablat*%CL_rad2deg), header= "Latitude (deg) - " + str(asc_desc)); 

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