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

// -----------------------------------------------------------
//> Time needed to access a location at Equator:  
//> Either ascending or descending passes are considered (but not both).  
//> The time is rounded to an integer number of days (days are rigorously 
//> numbers of Earth revolutions / ascending node). 
//> Plot parameter: 
//> 1: Maximum gap between ground tracks longitudes at equator
//> 2: View angle from the satellite ("roll angle")
//> 3: Corresponding "centre angle" (from Earth centre)
//> 4: Corresponding "incidence" at target location
//> (See: CL_gm_visiParams voir more details on 2,3,4) 
//> NB: Meaningful if the orbit's period is "sufficiently" less than 
//> the rotation period of the Earth.   
// ----------------------------------------------

// ---------------------------------------------------------------------------
// View angle from separation between ground tracks at equator
// => Should be a function in CelestLab in the future. 
// ---------------------------------------------------------------------------
function [ang] = View_angle(sma, inc, dlong, kind)
  angc = CL_op_equatorialSwath("dlon2cen",sma,inc,dlong); 
  ang = CL_gm_visiParams(sma,%CL_eqRad,"cen",angc,kind);
endfunction


// ---------------------------------------------------------------------------
// Maximum gap between longitudes after a certain number of orbits. 
// lgap: 1xN : longitude gap between 2 consecutive asc. node for each sma
// nborb: 1xN : max number of orbits for each semi major axis
// dlmax: 1xN : max gap between longitudes
// fracpso : 1*nbsat : satellites' arguments of latitude
// ---------------------------------------------------------------------------
function [dlmax] = Max_longitude_gap(lgap, nborb, pso_sat)
  Norb = 0:max(nborb); 

  fracpso = CL_rMod(pso_sat - pso_sat(1), -%pi, %pi) / (2*%pi); 

  // dl0: successives longitudes for each sma (NxNorb) 
  dl0 = lgap' * Norb; 
  dl1 = lgap' * ones(Norb); 

  Nb1 = ones(lgap') * Norb; // number of orbits (NxNorb)
  Nb2 = nborb' * ones(Norb); // max number of orbits (NxNorb)
  K = find(Nb1 > Nb2); 

  dl = [];
  for i = 1:length(fracpso)
    dli = CL_rMod(dl0 - dl1 * fracpso(i), 0, 2*%pi); 
    dli(K) = 2*%pi; 
    dl = [ dl, dli]; 
  end

  dl = gsort(dl, 'c', 'i'); // longitudes increasing order
  dl = [zeros(lgap)', dl, ones(lgap)' * 2*%pi]; // 0 and 2*%pi added (first longitude)

  dlmax = (max(dl(:,2:$) - dl(:,1:$-1), 'c'))';
endfunction


// ===========================================================================
// MAIN 
// ===========================================================================

h_min = 400.e3;  // min altitude
h_max = 1200.e3; // max altitude

Q_min = 1; // min number of "days"
Q_max = 10; // max number of "days"

sso = 1; // sun-synchronous
inc = 98*%CL_deg2rad; // inclination

ires = 2; // type of result

pso_sat = [0] * %CL_deg2rad; // in-plane satellite positions  

desc_param = list(..
   CL_defParam("Altitude of circular orbit - min", h_min, units=['m', 'km'], id='$hmin', valid='$hmin>=0' ),...
   CL_defParam("Altitude of circular orbit - max", h_max, units=['m', 'km'], id='$hmax', valid='$hmax>$hmin'  ),...
   CL_defParam("Maximum number of ''days'' (i.e. Earth revolutions)", Q_max, accv=1:20),...
   CL_defParam("Sun-synchronous orbit? (1=yes, 0=no)", sso, id='$sso', accv=[0,1]),...
   CL_defParam("Inclination (if not sun-synchronous)", inc, units=['rad', 'deg'], valid="$sso == 1 | ($x >0 & $x<180)"),...
   CL_defParam("Satellite(s) in-plane positions", pso_sat, units=['rad', 'deg'], dim=[1, 5]),...
   CL_defParam("Plot parameter (1=longitude gap, 2=sat angle, 3=centre angle, 4=incidence)", ires, accv=1:4)...
);

[h_min,h_max,Q_max,sso,inc,pso_sat,ires] = CL_inputParam(desc_param);

// ---------------------------------------------------------------------------

type_result = ["dl", "sat", "cen", "incid"];

nbsat = length(pso_sat);

nb = 100;
h = linspace(h_min, h_max, nb); // 100 altitude values

ecc = 0;

sma_min = h_min + %CL_eqRad;
sma_max = h_max + %CL_eqRad;

// add altitude corresponding to phased orbits
params = CL_op_searchRepeatOrbits(sma_min, sma_max, Q_min, nbsat*Q_max+1, ecc, sso, inc);
h1 = params(:,1)' - %CL_eqRad; // altitudes of "phased" orbits
h1 = h1(find(h1 > h_min & h1 < h_max)); // only altitudes in ]h_min, h_max[ ... in case

h = gsort([h,h1], 'c', 'i'); // all the altitudes (increasing order)

// inclination
sma = h + %CL_eqRad; // semi major axes
if (sso == 1)
  inc = CL_op_ssoJ2 ('i', sma, ecc); 
else
  inc = ones(sma) * inc;  // same size as sma
end

[lgap, T] = CL_op_paramsJ2(['lgap','nodper'],sma,ecc,inc); 
[drift_pom,drift_gom,drift_M]=CL_op_driftJ2(sma,ecc,inc); 
orbital_day = 2*%pi ./ (%CL_rotrBody - drift_gom); 
nborb_day = orbital_day ./ T; // number of orbits per orbital day

// results
res = zeros(Q_max, length(sma));   

for nbj = Q_min:Q_max
  nborb = round(nborb_day .* (nbj + 0.5)); // number of orbits during nbj+0.5 orbital days
  dlmax = Max_longitude_gap(lgap, nborb, pso_sat); 
 
  kind = type_result(ires); 
  if (kind == "dl")
    ang = dlmax;
  else
    ang = View_angle(sma,inc,dlmax/2,kind);
  end

  I = find(dlmax >= %pi);
  ang(I) = %pi/2; 

  res(nbj,:) = ang; 
end


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

f=scf();
f.visible = "on";
f.immediate_drawing = "off";
f.color_map = jetcolormap(max(3,Q_max+6));
f.color_map = f.color_map(4:$-3,:); 
//f.color_map = f.color_map(size(f.color_map,1):-1:1,:); // reverse order

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

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

for nbj=Q_max:-1:1

  plot2d([h, h($), h(1)]/1000, [res(nbj,:), %pi/2, %pi/2]  * %CL_rad2deg);
  CL_g_tag(a,nbj); 

  hdl = CL_g_select(a, "Polyline", nbj);
  hdl.closed = "on"; 
  hdl.fill_mode = "on"; 
  hdl.foreground = Gris; 
  hdl.background = nbj; 
end


y_titles = ["Max gap in longitude", "Satellite half view-angle", "View angle from Earth centre", "Incidence"] + " (deg)";

angmax = max(res);
a.data_bounds = [h_min/1000, 0; h_max/1000, min(angmax*1.05,%pi/2)*%CL_rad2deg];
a.tight_limits = "on"; 
a.title.text = "Accessibility at equator";
a.x_label.text = "Altitude (km)";
a.y_label.text = y_titles(ires);
CL_g_stdaxes();
CL_g_legend(a,string(Q_max:-1:Q_min) + " ", header="Access time (days)");


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