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

// -----------------------------------------------------------
//> Comparison between the DV of a continuous tangential thrust and the DV 
//> of the Hohmann transfer (in both cases to transform a circular orbit into 
//> another circular orbit). 
//> NB: The continuous thrust must be small enough for the final eccentricy
//> to be (close to) 0. The formula used does not guarantee that.   
//
// Author: A. Lamy
// -----------------------------------------------------------


h1min = 0.e3;   // m
h1max = 2000.e3;  // m

h2min = 20000.e3;  // m
h2max = 50000.e3;  // m

desc_param = list(..
   CL_defParam("Initial altitude - min",  h1min, units=['m', 'km'], id='$h1min', valid='$h1min>=0' ),..
   CL_defParam("Initial altitude - max", h1max, units=['m', 'km'], id='$h1max', valid='$h1max>$h1min'  ),..
   CL_defParam("Final altitude - min", h2min, units=['m', 'km'], id='$h2min', valid='$h2min>=0' ),..
   CL_defParam("Final altitude - max", h2max, units=['m', 'km'], id='$h2max', valid='$h2max>$h2min' )..
);

[h1min, h1max, h2min, h2max] = CL_inputParam(desc_param);

nbpts = 30;
tab_h1 = linspace(h1min,h1max,nbpts);
tab_h2 = linspace(h2min,h2max,nbpts+1);

ratio = [];
for h2 = tab_h2;
  dga_i = %CL_eqRad + tab_h1;
  dga_f = ones(tab_h1) * (%CL_eqRad + h2);

  // DV are >= 0  
  dv_hohmann = CL_man_dvHohmann(dga_i, dga_f); 
  dv_continuous = abs(sqrt(%CL_mu ./ dga_i) - sqrt(%CL_mu ./ dga_f)); 

  I = find(dv_hohmann == 0); 
  dv_hohmann(I) = %nan; // avoids division by 0
  
  r = 100 * (dv_continuous - dv_hohmann) ./ dv_hohmann; 
  
  I = find(isnan(r)); 
  r(I) = 0;  
  
  [ratio] = [ratio; r];
   
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(ratio); // grayplot min/max
zmax=max(ratio);

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

contour2d(tab_h1/1000, tab_h2/1000, ratio', levels, style=Coul1*ones(levels));
CL_g_tag(a,1);

contour2d(tab_h1/1000, tab_h2/1000, ratio', sublevels, style=Coul2*ones(sublevels));
CL_g_tag(a,2);

// general setting
CL_g_stdaxes(a);
a.data_bounds = [h1min,h2min;h1max,h2max]/1000;
a.tight_limits="on";
a.title.text = "Relative increase of |DV continuous| wrt |DV hohmann| (%)";
a.x_label.text = "Initial altitude (km)";
a.y_label.text = "Final 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(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";