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

// -----------------------------------------------------------
//> Cost of change of perigee/apogee altitude:
//> Total cost (DV >=0) needed for the transfer
//
//> The transfer consists in 2 maneuvers (no intermediate orbit). 
//> NB: The names "apogee" and "perigee" for the final orbit are not 
//>     (necessarily) what they usually mean. 
//>     They should be understood as altitude on one side of the line 
//>     of apsides, and altitude on the opposite side. 
//
// Author: A. Lamy
// -----------------------------------------------------------

// -----------------------------------------------------------
// data input
// -----------------------------------------------------------
k_init = 1; 
desc_param = list(..
   CL_defParam("Default parameters: 1=GTO 2=LEO", k_init, accv=[1,2])..
); 
k_init = CL_inputParam(desc_param);

if (k_init == 1)
  hpi = 200.e3;   // initial perigee altitude - m 
  hai = 36000.e3;  // initial apogee altitude - m
  hpfmin = 200.e3;  // final "perigee" altitude - m
  hpfmax = 10000.e3;  // final "perigee" altitude - m
  hafmin = 20000.e3;  // final "apogee" altitude - m
  hafmax = 100000.e3;  // final "apogee" altitude - m
else
  hpi = 500.e3;   // initial perigee altitude - m 
  hai = 500.e3;  // initial apogee altitude - m
  hpfmin = 200.e3;  // final "perigee" altitude - m
  hpfmax = 2000.e3;  // final "perigee" altitude - m
  hafmin = 200.e3;  // final "apogee" altitude - m
  hafmax = 2000.e3;  // final "apogee" altitude - m
end


desc_param = list(..
   CL_defParam("Initial perigee altitude",  hpi, units=['m', 'km'], id='$hpi', valid='$hpi>=0' ),..
   CL_defParam("Initial apogee altitude", hai, units=['m', 'km'], id='$hai', valid='$hai>=$hpi'  ),..
   CL_defParam("Final perigee altitude - min", hpfmin, units=['m', 'km'], id='$hpfmin', valid='$hpfmin>=0' ),..
   CL_defParam("Final perigee altitude - max", hpfmax, units=['m', 'km'], id='$hpfmax', valid='$hpfmax>$hpfmin' ),..
   CL_defParam("Final apogee altitude - min", hafmin, units=['m', 'km'], id='$hafmin', valid='$hafmin>=0' ),..
   CL_defParam("Final apogee altitude - max", hafmax, units=['m', 'km'], id='$hafmax', valid='$hafmax>$hafmin' )..
);

[hpi,hai,hpfmin,hpfmax,hafmin,hafmax] = CL_inputParam(desc_param);


// -----------------------------------------------------------
// computation
// -----------------------------------------------------------
N = 100;
tab_hpf = linspace(hpfmin,hpfmax,N);
tab_haf = linspace(hafmin,hafmax,N);


res_dv = [];
for haf = tab_haf;

  rpf = %CL_eqRad + min(tab_hpf, haf); // radius at (true) perigee 
  raf = %CL_eqRad + max(tab_hpf, haf); // radius at (true) apogee
  rai = %CL_eqRad + hai; 
  rpi = %CL_eqRad + hpi; 

  [ai,ei] = CL_op_rarp2ae(rai,rpi);
  [af,ef] = CL_op_rarp2ae(raf,rpf);

  // Check 2 values for 1st maneuver
  // The final perigee is on the same side as the initial perigee
  [dv1] = CL_man_dvHohmannG(ai,ei,af,ef,posman1=0,rotation=0); // dv > 0
  [dv2]  = CL_man_dvHohmannG(ai,ei,af,ef,posman1=1,rotation=0); // dv > 0

  [res_dv] = [res_dv; min(dv1, dv2)];

end

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

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

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

Coul1 = 5;
Coul2 = 10;
Coul3 = 2;

zmin=min(res_dv); // grayplot min/max
zmax=max(res_dv);

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

contour2d(tab_hpf/1000, tab_haf/1000, res_dv', levels, style=Coul1*ones(levels));
CL_g_tag(a,1);

contour2d(tab_hpf/1000, tab_haf/1000, res_dv', sublevels, style=Coul2*ones(sublevels));
CL_g_tag(a,2);

// general setting
CL_g_stdaxes(a);
a.data_bounds = [hpfmin,hafmin;hpfmax,hafmax]/1000;
a.tight_limits="on";
a.title.text = "DV needed for perigee/apogee altitude change (m/s)";
a.x_label.text = "Final perigee altitude (km)";
a.y_label.text = "Final apogee altitude (km)";

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

h = CL_g_select(a,"Text", 1);
CL_g_set(h, "text", string(abs(strtod(h.text))));

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";