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

// -----------------------------------------------------------
//> Swing-by = change of excess velocity vector caused by a body 
//> fly-by (or swing-by). 
//> 
//> Computed quantities: 
//> - Norm of velocity increment: Norm of difference between departure 
//> excess velocity vector and arrival excess velocity vector. 
//> - Turn angle: angle between  between departure 
//> excess velocity vector and arrival excess velocity vector. 
//> 
//> The results can be computed for the following planets: 
//> Venus, Earth, Mars, Jupiter
//
// Author: A. Lamy
// -----------------------------------------------------------

// Excess velocity (m/s)
vinfmin = 0.e3; 
vinfmax = 10.e3; 

// Periapsis Altitude (m)
hpmin = 0; 
hpmax = 10000.e3; 

// Body number
bodynum = 3; 

// Plot option 
opt = 1; 

desc_param = list(..
   CL_defParam("Hyperbolic excess velocity - min", vinfmin, units=['m/s', 'km/s'], id='$vinfmin', valid='$vinfmin >= 0'), ..
   CL_defParam("Hyperbolic excess velocity - max", vinfmax, units=['m/s', 'km/s'], id='$vinfmax', valid='$vinfmax > $vinfmin'), ..
   CL_defParam("Periapsis altitude - min", hpmin, units=['m', 'km'], id='$hpmin', valid='$hpmin >= 0'), ..
   CL_defParam("Periapsis altitude - max", hpmax, units=['m', 'km'], id='$hpmax', valid='$hpmax > $hpmin'), ..
   CL_defParam("Planet (2=Venus, 3=Earth, 4=Mars, 5=Jupiter)", bodynum, accv = 2:5), ..
   CL_defParam("Plot: 1=dv, 2=turn angle", opt, accv=1:2) ..
);

[vinfmin, vinfmax, hpmin, hpmax, bodynum, opt] = CL_inputParam(desc_param);

BodyNames = ["Venus", "Earth", "Mars", "Jupiter"]; 
BodyNums = [2, 3, 4, 5]; 

bodyname = BodyNames(find(bodynum == BodyNums)); 

// -----------------------
// Computation
// -----------------------

// planetary constants 
MU = CL_dataGet("body." + bodyname + ".mu"); 
R = CL_dataGet("body." + bodyname + ".eqRad");; 

// min value for vinf not 0 to avoid %nan 
vinfmin1 = max(vinfmin, vinfmax/10000); 
vinf = linspace(vinfmin1, vinfmax, 70);
hp = linspace(hpmin, hpmax, 70); 

[Vinf, Hp] = ndgrid(vinf, hp); 
all_vinf = matrix(Vinf, 1, -1); 
all_hp = matrix(Hp, 1, -1); 

[turnang, dv] = CL_ip_flybyParams("vinf", all_vinf, "rp", all_hp + R, ["turnang", "dv"], mu=MU); 

if (opt == 1)
  Z = matrix(dv / 1000, size(Vinf)); 
else
  Z = matrix(turnang * %CL_rad2deg, size(Vinf)); 
end
  
  
// -----------------------
// Plot
// -----------------------

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

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

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

contour2d(vinf/1.e3, hp/1.e6, Z, sublevels, style=color("lightblue")*ones(sublevels));
CL_g_tag(a, 2); 

contour2d(vinf/1.e3, hp/1.e6, Z, levels, style=color("blue") * ones(levels));
CL_g_tag(a, 1); 

// adjustments
CL_g_stdaxes(a);

a.data_bounds = [vinfmin/1.e3, hpmin/1.e6; vinfmax/1.e3, hpmax/1.e6];
a.tight_limits="on";

if (opt == 1)
  a.title.text = bodyname + " swing-by: Norm of velocity increment (km/s)";
else
  a.title.text = bodyname + " swing-by: Turn angle (deg)";
end

a.x_label.text = "Hyperbolic excess velocity (km/s)";
a.y_label.text = "Periapsis altitude (1.e3 km)";

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=color("blue");
h.font_size=3;
h.font_style=8;

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

f.immediate_drawing = "on";