//  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 Hohmann - Bi-elliptic (circular orbits):
//> Initial/final altitudes such that DV bi-elliptic < DV Hohmann.
//
// Author: A. Lamy
// -----------------------------------------------------------


// Default values
h1min = 0.e3;
h1max = 6000.e3;
h2min = 80000.e3;
h2max = 200000.e3;
hapot = 200000.e3;

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' ),..
   CL_defParam("Intermediate apogee altitude", hapot, units=['m', 'km'], id='$hapot', valid='$hapot>=0'  )..
);

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


nbpts = 80;
tab_h1 = linspace(h1min,h1max,nbpts);
tab_h2 = linspace(h2min,h2max,nbpts);

deltav = [];

for h2 = tab_h2;
  dga1 = %CL_eqRad + tab_h1;
  dga2 = (%CL_eqRad + h2) * ones(dga1);
  rapot   = (%CL_eqRad + hapot) * ones(dga1);

  dv_hohm = CL_man_dvHohmann(dga1, dga2);
  dv_biel = CL_man_dvBiElliptic(dga1, dga2, rapot);

  [deltav] = [deltav ; max(dv_hohm - dv_biel,0)];
end


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

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

f.color_map = jetcolormap(64);
Noir = addcolor([0,0,0]);
Gris = addcolor([1,1,1]*0.4);

a=gca();

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

colorbar(zmin,zmax,colminmax=[1,64],fmt="%.1f");
Sgrayplot(tab_h1/1.e6, tab_h2/1.e6, deltav',colminmax=[1,64],zminmax=[zmin,zmax]);

// general setting
a.data_bounds = [h1min,h2min;h1max,h2max]*1.e-6;
a.margins = [0.12,0.12,0.12,0.12];
CL_g_stdaxes();

a.title.text = "Gain of bi-elliptic transfer over Hohmann transfer (m/s)";
a.x_label.text = "Initial altitude (10^3 km)";
a.y_label.text = "Final altitude (10^3 km)";

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