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

// -----------------------------------------------------------
//> Secular drift due to J2 of a combination of variables: 
//> n1*dw/dt + n2*dW/dt + n3*dM/dt + n4*wE + n5*wS + n6*n 
//> with: 
//> n1..n6: Integers
//> dw/dt: Drift of argument of perigee
//> dW/dt: Drift of right ascension of ascending node
//> dM/dt: Drift of mean anomaly
//> wE:    Rotation rate of the Earth
//> wS:    (apparent) Rotation rate of the Sun
//> n:     Keplerian mean motion
//
// Author: A. Lamy
// -----------------------------------------------------------


smamin = 6400.e3;
smamax = 8000.e3;

incmin = 0 * %CL_deg2rad;
incmax = 180 * %CL_deg2rad;

ecc = 0;  

coef = [0, 0, 1, 0, 0, -1]; 

desc_param = list(..
   CL_defParam("Semi-major axis - min", smamin, units=['m', 'km'], id='$amin', valid='$amin>0'),..
   CL_defParam("Semi-major axis - max", smamax, units=['m', 'km'], id='$amax', valid='$amax>$amin'),..
   CL_defParam("Inclination - min", incmin, units=['rad', 'deg'], id='$imin', valid='$imin>=0 & $imin<=180'),..
   CL_defParam("Inclination - max", incmax, units=['rad', 'deg'], id='$imax', valid='$imax>$imin & $imax<=180'),..
   CL_defParam("Eccentricity", ecc, valid='$x>=0 & $x<1'),..
   CL_defParam("Coefficients (dw/dt, dW/dt, dM/dt, wE, wS, n)", coef, accv=-10:10, dim=6, valid='round($x) == $x')..
);

[smamin,smamax,incmin,incmax,ecc,coef] = CL_inputParam(desc_param);

nbpts = 80;
sma = linspace(smamin, smamax, nbpts);
inc = linspace(incmin, incmax, nbpts);
mm = CL_kp_params(['mm'], sma); // keplerian mean motion


res = []; // res(inc, sma) = raan
for i = inc;
  [drift_pom, drift_gom, drift_anm] = CL_op_driftJ2(sma,ecc,i);

  val = coef(1)*drift_pom + coef(2)*drift_gom + coef(3)*drift_anm + ..
        coef(4)*%CL_rotrBody + coef(5)*%CL_rotrBodySun + coef(6)*mm;   
               
  res = [res; val * %CL_rad2deg*86400]; // deg/j
end

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

// generates the name of the plotted parameter
function [str] = Name_param(coef)
  names = ["dw/dt", "dW/dt", "dM/dt", "wE", "wS", "n"]; 
  signes = [" -",""," +"]; 
  str = ""; 
  Inz = [find(coef <> 0), [7]]; 

  for k=Inz(1):6; 
    if (k > Inz(1) | coef(k)<0); str = str + signes(sign(coef(k))+2); end
    if (coef(k) <> 0)
      if (abs(coef(k)) <> 1); str = str + " " + sprintf("%d", abs(coef(k))) + " *" ; end
      str = str + " " + names(k); 
    end
  end
endfunction


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

f.color_map = jetcolormap(32);
Coul1 = 5;
Coul2 = 10;
Coul3 = 2;
Noir = addcolor([1,1,1]*0);

a=gca();

[levels, sublevels] = CL_autoLevels(min(res), max(res));

contour2d(sma/1000, inc*%CL_rad2deg, res', levels, style=Coul1*ones(levels));
CL_g_tag(a,1);

contour2d(sma/1000, inc*%CL_rad2deg, res', sublevels, style=Coul2*ones(sublevels));
CL_g_tag(a,2);

// general setting
CL_g_stdaxes(a);
a.data_bounds = [smamin/1000,incmin*%CL_rad2deg; smamax/1000,incmax*%CL_rad2deg];
a.tight_limits="on";
a.x_label.text = "Semi major axis (km)";
a.y_label.text = "Inclination (deg)";

a.title.text = Name_param(coef) + " (deg/day)" + "  -  ecc = " + string(ecc) ;

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

h = CL_g_select(a, "Text");
h.font_foreground=Coul3;
h.font_size=3;
h.font_style=8;

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

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