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

// -----------------------------------------------------------
//> Ellipse with radius and velocity vectors at equally 
//> distributed true anomalies. 
// 
//
// Author: A. Lamy
// -----------------------------------------------------------

// choice of input parameters
input_param = 2; 

desc_param = list(..
   CL_defParam("Input parameters: 1=sma/ecc, 2=ha/hp", input_param, accv=[1,2])..
);
input_param = CL_inputParam(desc_param);

if (input_param == 1)
  // sma/ecc
  sma = 24371.136e3; 
  ecc = 0.7300850; 
  desc_param = list(..
     CL_defParam("Semi major axis", sma,  units=['m', 'km'], valid="$x>0 & $x<=1.e6" ),..
     CL_defParam("Eccentricity",    ecc,  valid="$x>=0 & $x<1"  )..
  );
  [sma, ecc] = CL_inputParam(desc_param);
else
  // ha/hp
  ha = 35786.e3;
  hp = 200.e3;
  desc_param = list(..
     CL_defParam("Apogee altitude", ha, units=['m', 'km'], valid="$x>$hp & $x<1.e6" ),..
     CL_defParam("Perigee altitude", hp, units=['m', 'km'], id="$hp", valid="$x>=0" )..
  );
  [hp, ha] = CL_inputParam(desc_param);
  [sma, ecc] = CL_op_rarp2ae(%CL_eqRad + ha, %CL_eqRad + hp); 
end

 
anv1 = linspace(0, 2*%pi, 200);
anv2 = linspace(0, 2*%pi, 13);

kep1 = CL_kp_characteristics(sma,ecc,anv1);
kep2 = CL_kp_characteristics(sma,ecc,anv2);

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

a=gca();
a.title.text = "Semi-major axis = " + string(sma/1000) + " km, Eccentricity = " + string(ecc);
a.title.font_size=3;


// --------------------------------------
// plot Earth
// --------------------------------------

x1 = %CL_eqRad .* cos(anv1);
y1 = %CL_eqRad .* sin(anv1);

plot2d(x1,y1,style=2);

CL_g_tag(a,1);
h = CL_g_select(a, "Polyline",1);
h.background = 12;
h.fill_mode = "on";


// --------------------------------------
// plot ellipse
// --------------------------------------

r1 = kep1.r;
x1 = r1 .* cos(anv1);
y1 = r1 .* sin(anv1);

r2 = kep2.r;
x2 = r2 .* cos(anv2);
y2 = r2 .* sin(anv2);

plot2d(x1,y1,style=16);

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

// --------------------------------------
// plot velocity
// --------------------------------------

vel2 = (kep2.vel/kep2.vel(1)) * sma * 0.4;
gam2 = kep2.slope;

xp2 = x2 + vel2 .* cos(anv2+%pi/2-gam2);
yp2 = y2 + vel2 .* sin(anv2+%pi/2-gam2);

for i=1:length(anv2)-1;
  plot2d([x2(i), xp2(i)], [y2(i), yp2(i)], style=28);
end

CL_g_tag(a,4);
h = CL_g_select(a, "Polyline",4);
h.polyline_style=4;
h.thickness=2;
h.arrow_size_factor=0.7;


// --------------------------------------
// plot radius vectors
// --------------------------------------

for i=1:length(anv2)-1;
  plot2d([0, x2(i)], [0, y2(i)], style=16);
end

CL_g_tag(a,5);
h = CL_g_select(a, "Polyline",5);
h.polyline_style=4;
h.thickness=2;
h.arrow_size_factor=0.5;


a.isoview="on";
a.axes_visible=["off","off","off"];
a.margins = [0.1,0.1,0.12,0.1];
a.tight_limits = "on";
db = a.data_bounds;
db = db + [-sma, -sma; sma, sma] * 0.05;
a.data_bounds = db;

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


// --------------------------------------
// print information
// --------------------------------------

printf("\n %10s%-6s %10s%-6s %10s%-6s %10s%-6s\n", ...
       'True anom','(deg)', 'Radius', '(km)', 'Velocity', '(km/s)', 'Time', '(s)');

for i=1:length(anv2)
  printf(" %10.0f %18.1f %17.2f %15.1f\n", anv2(i)*180/%pi, ...
         kep2.r(i)/1000, kep2.vel(i)/1000, kep2.tperi(i));
end