//  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 chosen units are: 
//> - Body radius (R) for the periapsis radius, 
//> - Escape velocity at 1 body radius (=sqrt(2*MU/R)) for the excess velocity. 
//> Using these units makes the results body-independent. 
//> 
//> Approximate values for some planets: 
//>                     Radius (km)    Escape velocity (km/s)
//> Venus:                 6052            10.4
//> Earth:                 6378            11.2
//> Mars:                  3396             5.0
//> Jupiter:              71492            59.5
//> 
//> Note: 
//> It can be shown that the velocity increment at a distance r from 
//> the body center is maximum for an excess velocity equal to 
//> sqrt(mu/r) = escape velocity at the same position divided by sqrt(2). 
//> The maximum velocity increment is also equal to sqrt(mu/r). 
//> It follows that the maximum achievable velocity increment is 
//> sqrt(mu/R), with R = body radius. 
// 
// Author: A. Lamy
// -----------------------------------------------------------

// Excess velocity - unit = escape velocity
vinfmin = 0.; 
vinfmax = 1.; 

// Radius at periapsis - unit = body radius
rpmin = 1; 
rpmax = 3; 

// Plot option 
opt = 1; 


desc_param = list(..
   CL_defParam("Hyperbolic excess velocity - min [unit = escape velocity]", vinfmin, id='$vinfmin', valid='$vinfmin >= 0'), ..
   CL_defParam("Hyperbolic excess velocity - max [unit = escape velocity]", vinfmax, id='$vinfmax', valid='$vinfmax > $vinfmin'), ..
   CL_defParam("Periapsis radius - min [unit = body radii]", rpmin, id='$rpmin', valid='$rpmin >= 1'), ..
   CL_defParam("Periapsis radius - max [unit = body radii]", rpmax, id='$rpmax', valid='$rpmax > $rpmin'), ..
   CL_defParam("Plot: 1=dv, 2=turn angle", opt, accv=1:2) ..
);

[vinfmin, vinfmax, rpmin, rpmax, opt] = CL_inputParam(desc_param);


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

// Arbitrary choice of constants ! 
MU = 1; 
R = 1; 
ESCVEL = sqrt(2*MU/R); 

// min value for vinf not 0 to avoid %nan 
vinfmin1 = max(vinfmin, vinfmax/10000); 
vinf = linspace(vinfmin1, vinfmax, 80);
rp = linspace(rpmin, rpmax, 80); 

[Vinf, Rp] = ndgrid(vinf, rp); 
all_vinf = matrix(Vinf, 1, -1); 
all_rp = matrix(Rp, 1, -1); 

[turnang, dv] = CL_ip_flybyParams("vinf", all_vinf * ESCVEL, "rp", all_rp * R, ["turnang", "dv"], mu=MU); 

if (opt == 1)
  // Z = matrix(dv ./ (all_vinf * ESCVEL), size(Vinf)); 
  Z = matrix(dv / ESCVEL, 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, rp, Z, sublevels, style=color("lightblue")*ones(sublevels));
CL_g_tag(a, 2); 

contour2d(vinf, rp, Z, levels, style=color("blue") * ones(levels));
CL_g_tag(a, 1); 

// adjustments
CL_g_stdaxes(a);

a.data_bounds = [vinfmin, rpmin; vinfmax, rpmax];
a.tight_limits="on";

if (opt == 1)
  a.title.text = "Swing-by: Norm of velocity increment / escape velocity on body surface";
else
  a.title.text = "Swing-by: Turn angle (deg)";
end

a.x_label.text = "Hyperbolic excess velocity / escape velocity on body surface";
a.y_label.text = "Periapsis radius / body radius";

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