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

// -----------------------------------------------------------
//> Repeating ground tracks:  
//> Gives the semi-major axis of orbits such that the ground tracks 
//> repeat after a certain number of planet revolutions (Q). 
//
// Author: A. Lamy
// -----------------------------------------------------------

smamin = 6878.e3;
smamax = 7278.e3;
Qmin = 1;
Qmax = 20;
sso = 1; // sun-synchronous
inc = 0; // unused if SSO
ecc = 0; 
model = 2;
iplot = 1; // x axis = sma

desc_param = list(..
   CL_defParam("Semi-major axis - min", smamin, units=['m', 'km'], id='$smamin', valid='$smamin>=0' ),..
   CL_defParam("Semi-major axis - max", smamax, units=['m', 'km'], id='$smamax', valid='$smamax>$smamin' ),..
   CL_defParam("Eccentricity", ecc, valid='$x>=0 & $x <1' ),..
   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, units=['rad', 'deg'], valid="$sso == 1 | ($x >=0 & $x<=180)"),..
   CL_defParam("Q (number of planet revolutions / node) - min", Qmin, id='$Qmin', accv = 1:100),..
   CL_defParam("Q (number of planet revolutions / node) - max", Qmax, id='$Qmax', accv = 1:100, valid='$Qmax >= $Qmin'),..
   CL_defParam("Model (1=kepler, 2=J2)", model, accv = [1,2], id='$model'),..
   CL_defParam("X_axis (1=semi-major axis, 2=altitude)", iplot, accv = [1,2])..
);

[smamin,smamax,ecc,sso,inc,Qmin,Qmax,model,iplot] = CL_inputParam(desc_param);

j2 = %CL_j1jn(2);
if (model == 1); j2 = 0; end

result = CL_op_searchRepeatOrbits(smamin,smamax,Qmin,Qmax,ecc,sso,inc,j2=j2);


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

if (iplot == 1) 
  xaxis = result(:,1); 
  xlabel = "Semi major axis [km]";
  xaxis_min = smamin; 
  xaxis_max = smamax; 
else
  xaxis = result(:,1) - %CL_eqRad; 
  xlabel = "Altitude (semi major axis minus equ. radius) [km]";
  xaxis_min = smamin - %CL_eqRad; 
  xaxis_max = smamax - %CL_eqRad; 
end

plot(xaxis/1000, result(:,6),'d');  
[xmin, xmax] = CL_graduate(min(xaxis), max(xaxis));

a = gca();
a.data_bounds=[max(xmin,xaxis_min)/1000,Qmin-1;min(xmax,xaxis_max)/1000,Qmax+1];
a.tight_limits = "on";
CL_g_stdaxes(a);

a.grid_position = "background";
a.x_label.text = xlabel;
a.y_label.text = 'Number of planet revolutions (Q)';
if (sso == 0)
  a.title.text = 'Repeat orbits (Ecc = ' + string(ecc) + ', Inc = ' + string(inc*%CL_rad2deg) + ' deg)';
else
  a.title.text = 'Sun-synchronous repeat orbits (Ecc = ' + string(ecc) + ')';
end

h = CL_g_select(a, "Polyline");
h.mark_size=6;
h.mark_background=18;
h.mark_foreground=2;

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