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

// -----------------------------------------------------------
//> True anomaly as function of mean anomaly  
//> for elliptical or hyperbolic orbits. 

// Author: AL
// -----------------------------------------------------------


// Defaut values
M_min = 0; 
M_max = 2*%pi; 
tab_ecc = [0 0.2 0.5 0.7 0.9 0.99 1.01 1.1 1.5 2.5 5 20];
y_axis = 2; 

desc_param = list(..
   CL_defParam("Mean anomaly - min", M_min, units=['rad', 'deg'], id = '$M_min' ),...
   CL_defParam("Mean anomaly - max", M_max, units=['rad', 'deg'], id = '$M_max', valid='$M_max > $M_min' ),...
   CL_defParam("Eccentricity values", tab_ecc, dim=[1, 20], valid='$x >= 0 & abs($x-1) > 1.e-5'),...
   CL_defParam("y-axis (1=true anomaly, 2=true minus mean)", y_axis, accv=[1,2])...
);

[M_min,M_max,tab_ecc,y_axis] = CL_inputParam(desc_param);


nbpts = 1200;
tab_M = linspace(M_min,M_max,nbpts);


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

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

nb = length(tab_ecc);
f.color_map = jetcolormap(min(max(nb,3),100));

Noir = addcolor([0,0,0]);
GrisF = addcolor([1,1,1]*0.4);

ymin = +%inf; 
ymax = -%inf;

for k = 1:length(tab_ecc);
  ecc = tab_ecc(k);
  tab_v = CL_kp_M2v(ecc, tab_M);

  if (y_axis == 1) 
     y = tab_v;
  else
     y = tab_v - tab_M; 
  end

  plot2d(tab_M*%CL_rad2deg, y*%CL_rad2deg, style=k);
  ymin = min(min(y), ymin); 
  ymax = max(max(y), ymax); 

end

a=gca();

CL_g_legend(a, string(tab_ecc), header="Eccentricity");

// adjustments

// graduations
[xmin, xmax, nb] = CL_graduate(M_min*%CL_rad2deg, M_max*%CL_rad2deg); 
xticks = linspace(xmin,xmax,nb+1);
a.x_ticks = tlist(['ticks', 'locations', 'labels'], xticks, string(xticks));

[ymin, ymax] = CL_graduate(ymin*%CL_rad2deg,ymax*%CL_rad2deg); 
a.data_bounds = [xmin,ymin; xmax,ymax];
a.tight_limits = "on"; 

CL_g_stdaxes(a, colg=GrisF);

a.x_label.text = "Mean anomaly (deg)";
if (y_axis == 1)
  a.title.text = "True anomaly";
  a.y_label.text = "Angle (deg)";
else
  a.title.text = "True anomaly minus mean anomaly";
  a.y_label.text = "Angle (deg)";
end


// proprietes des courbes
h=CL_g_select(a, "Polyline");
h.thickness=2;

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