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

// ----------------------------------------------------------- 
//> Evaluation of space/time coverage for a constellation of satellite: 
//> passes of spacecraft in a constellation at a given latitude. 
//> 
//> Each rectangle in the graph corresponds to a particular pass: 
//> X-axis : longitude of pass, Y-axis: time of pass. 
//> The 2 input parameters "reference size" and "reference time length" define
//> the size of the rectangle along the X and Y axes respectively. 
//> 
//> The orbits are circular, share the same repeat orbit parameters (N,P,Q), and are 
//> characterized by particular values of PSO (argument of latitude) and RAAN at initial 
//> time: 
//> - PSO index = k => pso = k * 2*pi/Q 
//> - RAAN index = l => raan = l * 2*pi/(N*Q+P)  
//> If k and l are integers, the ground track is the same as for the reference satellite
//> (defined by pso = 0 and raan = 0). 
//> 
// Author : A. Lamy
// -----------------------------------------------------------

// Notes
// By definition: 
// n / (we-gomdot) = N + P/Q = Norbcycle/Q; 

// delta_pso to have same tracks: 
// dlon = - (dpso / n) * (we - gomdot) = - dpso * Q / Norbcycle = K * 2*%pi / Norbcycle; 
// => dpso = - K * (2*%pi) / Q 
// => dlon =  (K * (2*%pi) / Q) * Q / Norbcycle = K * (2*%pi) / Norbcycle


// ----------------------- 
// Inputs
// ----------------------- 

N = 12; 
P = 7; 
Q = 10; 
NPQ = [N, P, Q]; 
sso = 0; 
inc = 66 * %CL_deg2rad; 
ecc = 0; 

// 1 = asc, -1 = desc, 0 = both
asc_desc = 1; 
lat = 0 * %CL_deg2rad; 

dx = 100.e3; 
dt = 2 * 86400; 

kpso_sat = [0, 5];   //  * 2*%pi/Q 
kgom_sat = [0, 0];  // * 2*%pi/Norbcycle

imodel = 2; // J2

// npq = [N, P, Q] 
function [ok] = valid_npq(NPQ)
  [N,P,Q] = (NPQ(1), NPQ(2), NPQ(3)); 
  ok = ( (N == round(N) & N > 0) & ..
         (P == round(P) & P >= 0 & P < Q) & ..
         (Q == round(Q) & Q > 0 & ((P == 0 & Q == 1) | CL_gcd(P, Q) == 1)) ); 
endfunction  

function [ok] = valid_lat(NPQ,ecc,sso,inc,lat)
  [sma, inc] = CL_op_repeat2smaInc(NPQ(1),NPQ(2),NPQ(3),ecc,sso,inc); 
  // adjust max latitude. 
  latmax = min(inc, %pi-inc) - 1.e-5;
  ok = (lat >= 0 & lat <= latmax); 
endfunction  

desc_param = list(..
   CL_defParam("Repeat orbit parameters (N,P,Q)", NPQ, dim=3, id='$NPQ', valid='valid_npq($x)' ),..
   CL_defParam("Sun-synchronous orbit? (1=yes, 0=no)", sso, id='$sso', accv=[0,1], valid="$sso==0 | $model==2"),..
   CL_defParam("Inclination (if not Sun-synchronous)", inc, id="$inc", units=['rad', 'deg'], valid="$sso == 1 | ($x >=0 & $x<=180)"),..
   CL_defParam("Latitude", lat, id='$lat', units=["rad", "deg"], valid="valid_lat($NPQ,0,$sso,$inc*%pi/180,$lat*%pi/180)"), ..
   CL_defParam("Which ground tracks? (1=ascending, -1=descending, 0=both)", asc_desc, accv = -1:1),..
   CL_defParam("Reference size on the ground along parallels", dx, units=["m", "km"], valid="dx > 0"),..
   CL_defParam("Reference time length", dt, units=["s", "day"], valid="dt > 0"),..
   CL_defParam("PSO indices (multiple of 2*pi/Q) for all spacecraft", kpso_sat, dim=[1, 20]),..
   CL_defParam("RAAN indices (multiple of 2*pi/(N*Q+P)) for all spacecraft", kgom_sat, dim=[1, 20]),..
   CL_defParam("Model (1=kepler, 2=J2)", imodel, id="$model", accv = [1,2])..
);

[NPQ,sso,inc,lat,asc_desc,dx,dt,kpso_sat,kgom_sat,imodel] = CL_inputParam(desc_param);
[N,P,Q] = (NPQ(1), NPQ(2), NPQ(3)); 

// J2 <- 0 if Keplerian model
if (imodel == 1)
  %CL_j1jn = %CL_j1jn; 
  %CL_j1jn(2) = 0; 
end

// Force same dimensions for kpso_sat and kgom_sat + removes unnecessary zeros
n = max(length(kpso_sat), length(kgom_sat))+1; 
kpso_sat(1,$+1:n) = 0; 
kgom_sat(1,$+1:n) = 0; 

while (kpso_sat($) == 0 & kgom_sat($) == 0 & length(kpso_sat) > 1)
  kpso_sat($) = []; 
  kgom_sat($) = []; 
end

nbsat = length(kpso_sat); 


// ----------------------- 
// Useful functions
// ----------------------- 

function [data] = compute_data_orb(N,P,Q,ecc,sso,inc); 
  [data.sma, data.inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso,inc); 
  // time derivatives
  [data.pomdot, data.gomdot, data.anmdot] = CL_op_driftJ2(data.sma,ecc,data.inc); 
  // mean motion$
  data.n = (data.pomdot + data.anmdot); 
  // period
  data.T = 2*%pi / data.n; 
  
  // Number of orbits per repeat, duration of cycle (sec)
  data.Norbcycle = N*Q+P; 
  data.Tcycle = data.T * data.Norbcycle; 

  // delta longitude after 1 orbit
  data.DL = -2*%pi / (N + P/Q); 
  
  // delta longitude in mesh 
  data.dLcycle = 2*%pi / data.Norbcycle; 
endfunction


// longitude at ascending node = 0
// asc => ascending, desc => descending
// norbd => N+P/Q
function [l_asc, t_asc, l_desc, t_desc] = compute_dlon(inc, lat, n, norbd); 
  omega = n / norbd; // we - gomdot

  // 1st longitude at latitude = lat (ascending) - pso in [-pi, pi]
  [pso, rra]  = CL_gm_orbSphTriangle(inc, "decl_a", lat, output=["aol", "rra"]);
  t_asc = pso / n; 
  l_asc = rra - t_asc * omega; 

  // 1st longitude at latitude = lat (descending) - pso in [-pi, pi]
  [pso, rra]  = CL_gm_orbSphTriangle(inc, "decl_d", lat, output=["aol", "rra"]);
  t_desc = pso / n; 
  l_desc = rra - t_desc * omega; 
endfunction

// xmin; xmax; ymin; ymax; color; t (rad / sec)
function [X] = gen_rects(l,t,dl,dt,icol)
  xmin = CL_rMod(l, -%pi, %pi) - dl/2; 
  xmax = CL_rMod(l, -%pi, %pi) + dl/2; 
  ymin = t - dt/2; 
  ymax = t + dt/2; 
  col = icol .* ones(l); 
  X = [xmin; xmax; ymin; ymax; col; t]; 
endfunction

function plot_rects(X, tabcol)
  // filter unvisible: 
  // I = find(X(1,:) > abs(DL/2) | X(2,:) < -abs(DL/2) | X(3,:) > Tcycle | X(4,:) < 0); 
  // X(:,I) = []; 

  X = CL_sortMat(X, X(6,:)); // date
  ux = %CL_rad2deg; 
  uy = 1 ./ 86400; 
  // x(top, left); y(top, left); w, h
  xrects([X(1,:)*ux; X(4,:)*uy; (X(2,:)-X(1,:))*ux; (X(4,:)-X(3,:))*uy], tabcol(X(5,:))); 
endfunction

function plot_gridx(x, ymin, ymax, icol)
  Y = [ymin, ymax]; 
  X = x' * ones(Y); 
  plot(X, Y); 
  h = CL_g_select(gce(), "Polyline"); 
  h.foreground = icol;
endfunction


function plot_gridy(xmin, xmax, y, icol)
  X = [xmin, xmax]; 
  Y = y' * ones(X); 
  plot(X, Y, "-"); 
  h = CL_g_select(gce(), "Polyline"); 
  h.foreground = icol;
  h.line_style = 8; 
endfunction


// ----------------------- 
// Computation
// ----------------------- 
data = compute_data_orb(N,P,Q,ecc,sso,inc); 

kpso_sat = CL_rMod(kpso_sat, -Q/2, Q/2); 
kgom_sat = CL_rMod(kgom_sat, -data.Norbcycle/2, data.Norbcycle/2); 

dt_sat = -kpso_sat * (2*%pi/Q) / data.n; 
dl_sat = kgom_sat * (2*%pi) / (N*Q+P) + kpso_sat * (2*%pi) / (N*Q+P); 

// dimension in "x"
dl = dx / (%CL_eqRad * cos(lat)); 

// gap in longitude at latitude lat / latitude 0
[dlon_asc, dt_asc, dlon_desc, dt_desc] = compute_dlon(data.inc, lat, data.n, N+P/Q);

// longitudes of ascending/descending tracks at latitude lat
// 3 repeat cycles
iorb = -data.Norbcycle : 2*data.Norbcycle-1; 

l_asc = CL_rMod(dlon_asc + iorb * data.DL, -%pi, %pi);
t_asc = dt_asc + iorb * data.T; 

l_desc = CL_rMod(dlon_desc + iorb * data.DL, -%pi, %pi);
t_desc = dt_desc + iorb * data.T; 

X = []; 
for k = 1:nbsat
  // asc
  if (asc_desc >= 0)
    icol = k; 
    X = [ X, gen_rects(l_asc+dl_sat(k), t_asc+dt_sat(k), dl, dt, icol) ]; 
 end
  // desc
  if (asc_desc <= 0)
    icol = k + nbsat; 
    X = [ X, gen_rects(l_desc+dl_sat(k), t_desc+dt_sat(k), dl, dt, icol) ]; 
  end
end

// filter rectangles (visible only)
I = find(X(1,:) > abs(data.DL/2) | X(2,:) < -abs(data.DL/2) | X(3,:) > data.Tcycle | X(4,:) < 0); 
X(:,I) = []; 

bounds = [-abs(data.DL/2), abs(data.DL/2), 0, data.Tcycle]; 


// ------------------
// Plot
// ------------------
f = scf(); 
f.visible = "on"; 
f.immediate_drawing = "off"; 
clf(); 

plot(0,0); 

// color indices
color_base = [2, 13, 22, 19, 27, 12, 28, 32, 17, 31]; 
color_sat_asc = repmat(color_base, 1, ceil(nbsat/10)); 
color_sat_desc = addcolor(0.5 + f.color_map(color_sat_asc,:) * 0.5); 

plot_rects(X, [color_sat_asc(1:nbsat), color_sat_desc(1:nbsat)]); 

plot(0,0); 
plot_gridy(-(data.DL/2)*%CL_rad2deg,(data.DL/2)*%CL_rad2deg, 0 : floor(data.Tcycle/86400), color("grey30")); 
plot_gridx(((-Q-1:Q+1) * data.dLcycle + CL_rMod(dlon_asc,-data.dLcycle/2,data.dLcycle/2)) * %CL_rad2deg, 0, data.Tcycle/86400, color_sat_asc(1)); 

a = gca(); 
CL_g_stdaxes(a); 
a.data_bounds = [-abs(data.DL/2)*%CL_rad2deg, 0; abs(data.DL/2)*%CL_rad2deg, (data.Tcycle)/86400]; 
a.tight_limits = "on"; 
a.grid = [-1,-1]; 

a.title.text = msprintf("N,P,Q = (%d,%d,%d), Inc = %.1f deg, Lat = %.1f deg", ..
               N, P, Q, data.inc*%CL_rad2deg, lat*%CL_rad2deg); 
a.x_label.text = "Longitude (deg)"; 
a.y_label.text = "Time (days)"; 

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