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

// -----------------------------------------------------------
//> DV needed to change the (elliptical) initial orbit into an orbit with 
//> a low enough (adjustable) perigee. 
//> Only one manoeuvre (at the apogee) is performed . 
//
// Author: A. Lamy
// -----------------------------------------------------------
 
ha =  800.e3; 

hpimin = 400.e3;
hpimax = 900.e3;

hpmin = 300.e3;
hpmax = 650.e3;

desc_param = list(..
   CL_defParam("Apogee altitude", ha, units=['m', 'km'], id='$ha', valid='$ha>=0' ),..
   CL_defParam("Perigee altitude of initial orbit - min", hpimin, units=['m', 'km'], id='$hpimin', valid='$hpimin>=0' ),..
   CL_defParam("Perigee altitude of initial orbit - max", hpimax, units=['m', 'km'], id='$hpimax', valid='$hpimax>$hpimin' ),..
   CL_defParam("Perigee altitude of final orbit - min", hpmin, units=['m', 'km'], id='$hpmin', valid='$hpmin>=0' ),..
   CL_defParam("Perigee altitude  of final orbit - max", hpmax, units=['m', 'km'], id='$hpmax', valid='$hpmax>$hpmin' )..
);

[ha,hpimin,hpimax,hpmin,hpmax] = CL_inputParam(desc_param);


nbpts = 40; 
tab_hpi = linspace(hpimin, hpimax, nbpts);
tab_hp = linspace(hpmin, hpmax, nbpts);


// -----------------------------------------------------------
// results
// -----------------------------------------------------------

tab_deltav = [];

for k = 1:length(tab_hp);

  // initial orbit
  [ai, ei] = CL_op_rarp2ae(%CL_eqRad + ha, %CL_eqRad + tab_hpi); 
  af = CL_op_rarp2ae(%CL_eqRad + ha, %CL_eqRad + tab_hp(k)); 

  [_X_, dv] = CL_man_dvSma(ai, ei, af, posman = "apo"); // returned norm ignored
  deltav = -dv(2,:); // signed: >0 if sma decreases
  
  tab_deltav = [tab_deltav; deltav];

end

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

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

Nmap = 256;
f.color_map = 0.95 + 0.05*jetcolormap(Nmap);

nb = length(tab_hp);
Noir = addcolor([0,0,0]);
GrisF = addcolor([1,1,1]*0.3);
GrisC = addcolor([1,1,1]*0.7);

a=gca();
CL_g_tag(a, 0);

Sgrayplot(tab_hpi/1000, tab_hp/1000, tab_deltav',...
          colminmax=[round(Nmap*0.2),round(Nmap*0.9)],...
          zminmax=[min(max(tab_deltav,0)),max(tab_deltav)], colout=[0,0]);

[levels, sublevels] = CL_autoLevels(min(max(tab_deltav,0)), max(tab_deltav), 8, 5);

contour2d(tab_hpi/1000, tab_hp/1000, tab_deltav', sublevels, style=GrisC*ones(sublevels));
CL_g_tag(a,2);

contour2d(tab_hpi/1000, tab_hp/1000, tab_deltav', levels, style=GrisF*ones(levels));
CL_g_tag(a,1);

// general setting
CL_g_stdaxes(a)
a.data_bounds = [hpimin, min(tab_hp); hpimax, max(tab_hp)]/1000;
a.tight_limits = "on";
a.title.text = "Deorbit cost (m/s) - Apogee altitude = " + sprintf("%.1f", ha/1000) + " km";
a.x_label.text = "Perigee altitude of initial orbit (km)";
a.y_label.text = "Perigee altitude of final orbit (km)";

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

h = CL_g_select(a, "Text");
CL_g_set(h, "text", string(strtod(h.text)));
h.font_foreground=Noir;
h.font_size=3;
h.font_style=8;

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

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