//  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), length(mltan)); 
  t = t .* ones(1,N); 
  x = x .* ones(1,N); 
  mltan = mltan .* 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; 
hlmin = 6; 
hlmax = 18;
lat = 40 * %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 - min", hlmin, units=['h'], id='$hlmin', valid='$hlmin>=0 & $hlmin<=24'),..
   CL_defParam("Mean local time of ascending node - max", hlmax, units=['h'], id='$hlmax', valid='$hlmax>=0 & $hlmax<=24 & $hlmax > $hlmin'),..
   CL_defParam("Latitude", lat, units=['rad', 'deg'], valid='abs($x) <= 90 & abs($x) <= min($inc, 180-$inc)'),..
   CL_defParam("Ground tracks (1=ascending, 2=descending)", asc_desc, accv=[1,2])..
);


[inc, hlmin, hlmax, lat, 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; 

tabhl = linspace(hlmin, hlmax, 60); 
tabs = linspace(0,1,60); // time in year (between 0 and 1)
tabt = t0 + tabs * yearlen; 

[Tabt, Tabhl] = ndgrid(tabt, tabhl); 
res = sat_info(inc, matrix(Tabhl,1,-1), lat, asc_desc, matrix(Tabt,1,-1)); 
sza = matrix(res.sza, size(Tabt));


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

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

nbcol = 64;
f.color_map = 0.5 + 0.5*jetcolormap(nbcol);

Noir = addcolor([0,0,0]); 
Gris = addcolor([1,1,1] * 0.3); 

a=gca();
CL_g_tag(a, 0);

x = tabs; 
y = tabhl; 
z = sza * %CL_rad2deg; 
zmin = min(z); 
zmax = max(z); 

colorbar(zmin,zmax,colminmax=[1,64],fmt="%.1f");

Sgrayplot(x, y, z, colminmax=[1,nbcol], zminmax=[zmin, zmax]);

[levels, sublevels] = CL_autoLevels(zmin, zmax);

contour2d(x, y, z, sublevels, style=Gris*ones(sublevels));
CL_g_tag(a,2);

contour2d(x, y, z, levels, style=Noir*ones(levels));
CL_g_tag(a,1);


// plot adjustments

h = CL_g_select(a, "Text", 2);
CL_g_delete(h);

h = CL_g_select(a, "Text");
CL_g_set(h, "text", string(strtod(h.text)));
h.font_foreground=Noir;
h.font_size=2;
h.font_style=8;

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

// titles
a.x_label.text = "Day / Month (" + string(year) + ")"; 
a.y_label.text = "Mean local time of ascending node (deg)"; 

hl_cal = CL_dat_cjd2cal(tabhl/24); 
str = ["asc", "desc"]; 
a.title.text = sprintf("Solar zenith angle (deg) - inc = %.1f deg, lat = %.1f deg, %s", inc*%CL_rad2deg, lat*%CL_rad2deg, str(asc_desc)); 
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,1], [labels, "1/1"]);
 
a.data_bounds = [0, hlmin; 1, hlmax]; 
a.tight_limits = "on"; 

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


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