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

// -----------------------------------------------------------
//> Orbit maintenance (compensation of atmospheric drag)  
//> 
//> The semi major-axis is supposed to decrease at a constant rate. 
//> Periodic semi major-axis corrections are then performed in 
//> order to keep the ground track or the argument of latitude inside
//> a given "window". 
//> The window size (w) is defined as a size on the ground for the longitudes
//> and a size at the altitude of the satellite the the argument of latitude. 
//> By definition, the controlled parameters must stay in [-w, w]. 
//> NB: 
//> - The orbit is circular. 
//> - The result is very theoretical: no uncertainties (on prediction of 
//> atmospheric drag, maneuver execution, ...) are taken into account. 
//
// Author: A. Lamy
// -----------------------------------------------------------

alt = 700.e3;
inc = 98 * %CL_deg2rad; 
ecc = 0; // fixed
absdadtmin = 1 / 86400; 
absdadtmax = 10 / 86400; 

sktype = 2; // longitudes
wsize = [1, 5, 10] *  1000; // window size on the ground (m)

desc_param = list(..
   CL_defParam("Altitude (semi-major axis minus equ. radius)", alt, units=['m', 'km'], valid='$x>=0 & $x < 5.e4'),..
   CL_defParam("Inclination", inc, units=['rad', 'deg'], valid="$x >=0 & $x<=180"),..
   CL_defParam("Sma decrease rate - min", absdadtmin, id="$absdadtmin", units=['m/s', 'm/day'], valid="$x > 0"),..
   CL_defParam("Sma decrease rate - max", absdadtmax,  units=['m/s', 'm/day'], valid="$x > $absdadtmin"), .. 
   CL_defParam("Type of station keeping: 1=arg of latitude, 2=ground track", sktype, accv=1:2), ..
   CL_defParam("Window size (on Earth surface)", wsize, units=['m','km'], valid='$x >= 0', dim=[1, 20]) ..
);

[alt, inc, absdadtmin, absdadtmax, sktype, wsize] = CL_inputParam(desc_param);

sma = %CL_eqRad + alt;
wsize = matrix(wsize, -1, 1); // column

// -----------------------------------------------------------
// computation
// -----------------------------------------------------------

absdadt = linspace(absdadtmin, absdadtmax, 100); 
dadt = -absdadt; 

// d(pso)/d(sma
[pomdot,gomdot,anmdot,dpomdotdaei,dgomdotdaei,danmdotdaei] = CL_op_driftJ2(sma,ecc,inc); 

psodot = pomdot + anmdot; 
dpsodotda = dpomdotdaei(1) + danmdotdaei(1); 

// A = deg 2 coefficient of pso model (rad/s^2) 
// A ~= -3/4 * (n/a) * da/dt
A = 0.5 * dpsodotda * dadt;

if (sktype == 1)
  wsize_rad = wsize / sma; // window size at satellite altitude (rad)
  K = A; 
else
  wsize_rad = wsize / %CL_eqRad; // window size on the ground (rad)
  // K ~ A * (omega_Earth / n)
  K = ((%CL_rotrBody - gomdot) / psodot) * A; 
end

// Time between maneuvers (sec)
Tman = sqrt(8 * wsize_rad * (1 ./ abs(K))); 


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

nbcolors = max(5, size(wsize, "*")); 
f.color_map = jetcolormap(nbcolors); 

for k = 1 : size(wsize,"*")

  plot(abs(dadt) * 86400, Tman(k,:) / 86400); 
  h = CL_g_select(gce(), "Polyline"); 
  h.thickness = 2; 
  h.foreground = k; 

end

a = gca();

a.x_label.text = "Semi major-axis decrease rate [m/day]";
a.y_label.text = "Time beween maneuvers [days]";

if (sktype == 1)
  a.title.text = "Station keeping (atmos. drag) - In-track position on the ground";
else
  a.title.text = "Station keeping (atmos. drag) - Ground track at equator";
end

CL_g_legend(a, string(wsize/1000), header="Window size (km)"); 

a.data_bounds = [abs(dadt(1))*86400, a.y_ticks.locations(1); abs(dadt($))*86400, a.y_ticks.locations($)]; 
a.tight_limits = "on"; 

CL_g_stdaxes(a);

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