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

// -----------------------------------------------------------
//> Eccentricity of frozen orbit:
//> Is such that the mean values of eccentricity and 
//> argument of perigee remain constant over time.
//> The quantity shown is ey = e*sin(w), w = argument of perigee. 
//> ey > 0 corresponds to: w = +90 degrees. 
//> ey < 0 corresponds to: w = -90 degrees. 
//>
//> NB: The relative accuracy decreases as the number of 
//> zonal terms increases because of numerical errors. 
//
// Author: A. Lamy
// -----------------------------------------------------------


smamin = 6800.e3;
smamax = 7800.e3;
incmin = 90 * %CL_deg2rad;
incmax = 110 * %CL_deg2rad;
nbz = 20;
nbzmax = length(%CL_j1jn); // max number of zonal coefs

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("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("Number of zonal coefs", nbz, accv=3:nbzmax)..
);

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

nbpts = 50;
sma = linspace(smamin, smamax, nbpts);
inc = linspace(incmin, incmax, nbpts+1); // +1 : used to matrices sizes


res = [];
for i = inc;
  [ecc,pom] = CL_op_frozenOrbit(sma,i*ones(sma),j1jn=%CL_j1jn(1:nbz));
  res = [res; ecc.*sin(pom)*1000]; // ecc * sin(pom) = ey multiplied by 1000
end

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

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 = "Eccentricity of frozen orbit (e*sin(w) x 1.e3)";

// change some properties
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";