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

function [varargout] = CL_ip_flybyParams(type_par1, par1, type_par2, par2, output, mu)
// Conversion between various fly-by parameters (hyperbolic orbits)
//
// Calling Sequence
// [par_out1, par_out2, ...] = CL_ip_flybyParams(type_par1, par1, type_par2, par2, output [, mu])
// 
// Description
// <itemizedlist><listitem>
// <p>Given two fly-by parameters, the function computes other ones. </p> 
// <p>The available fly-by parameters are:</p>
// <p>- <b>vinf</b>: hyperbolic excess velocity</p>
// <p>- <b>rp</b>: Radius at periapsis</p>
// <p>- <b>dinf</b>: Impact parameter = distance from the body center to the asymptote 
// of the hyperbola</p>
// <p>- <b>turnang</b>: Fly-by turn-angle = angle between arrival departure 
// excess velocity vectors</p>
// <p>- <b>vp</b>: Velocity at periapsis</p>
// <p>- <b>dv</b>: Norm of velocity increment</p>
// <p></p>
// <p><imageobject><imagedata fileref="flyby_params.png"/></imageobject></p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Note</b>:</p>
// <p>In the current version:</p>
// <p>- The 1st input parameter must be "vinf"</p>
// <p>- The 2nd input parameter must be "rp", "dinf" or "turnang"</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// type_par1: (string) Type of 1st input parameter. Only "vinf" accepted. (1x1)
// par1: Value of 1st input parameter [m, m/s, rad]. (1xN or 1x1)
// type_par2: (string) Type of 2nd input parameter: "rp", "dinf" or "turnang". (1x1)
// par2: Value of 2nd input parameter [m, m/s, rad]. (1xN or 1x1)
// output: (string) Names of the desired output parameters: "vinf", "rp", "dinf", "turnang", "vp" or "dv". (1xP)
// mu: (optional) Gravitational constant [m^3/s^2]. Default value is %CL_mu.
// par_out1, par_out2, ...: Values of the output parameters [m, m/s, rad].  (1xN) 
// 
// Authors
// CNES - DCT/SB
//
// See also
// CL_ip_flybyVectors
// 
// Bibliography
// H.D. Curtis, Orbital Mechanics for Engineering Students, ch. 8.
//
// Examples
// // Periapsis radius (m) and hyperbolic excess velocity (m/s)
// rp = 8.e6; 
// vinf = 3.e3; 
//
// // Flyby turn-angle (default value for "mu")
// turnang = CL_ip_flybyParams("vinf", vinf, "rp", rp, "turnang")
//
// // Check results
// CL_ip_flybyParams("vinf", vinf, "turnang", turnang, "rp") - rp  // => 0 
// 

// Declarations:

// ---------------------------
// Norm of velocity at infinity and impact parameter to periapsis radius.
// dinf = 0 => => rp = %nan
// NB: vinf must be <> 0
// ---------------------------
function [rp] = flybyPar_vinfd2rp(vinf, dinf, mu)
  I = find(dinf <= 0);
  dinf(I) = %nan;  
  rp = -mu ./ vinf.^2 + sqrt(((mu ./ vinf.^2).^2 + dinf.^2));
endfunction

// ---------------------------
// Norm of velocity at infinity and flyby turn angle angle to periapsis radius.
// turnang = 0 or %pi => dinf = %nan, rp = %nan 
// NB: vinf must be <> 0
// ---------------------------
function [rp] = flybyPar_vinfturnang2rp(vinf, turnang, mu)
  sintheta2 = sin(turnang / 2); 
  I = find(sintheta2 <= 0 | sintheta2 >= 1);
  sintheta2(I)= %nan;
  rp = (mu ./ vinf^2) .* (1 ./ sintheta2 - ones(turnang));
endfunction

// ---------------------------
// Norm of velocity at infinity and periapsis radius to various parameters
// rp == 0 => %nan
// NB: vinf must be <> 0
// ---------------------------
function [dinf, turnang, vp, dv] = flybyPar_vinfrp2all(vinf, rp, mu)
  I = find(rp <= 0);
  rp(I) = %nan;  
  
  dinf  = rp .* sqrt(1 + 2 * mu ./ (rp .* vinf.^2));
  sin_turnang_over_2 = 1 ./ (1 + rp .* (vinf.^2) / mu); 
  turnang = 2 * asin(sin_turnang_over_2); 
  vp = (vinf .* dinf) ./ rp; 
  dv = 2 * vinf .* sin_turnang_over_2; 
endfunction


// Code

if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// ------------------
// Argument Checking
// ------------------
// number of argument (left/right)
[lhs, rhs] = argn();

if (rhs < 4)
  CL__error("Invalid number of input arguments (at least 4 expected)");
end

if (typeof(type_par1) <> "string" | size(type_par1, "*") <> 1 | ..
    typeof(type_par2) <> "string" | size(type_par2, "*") <> 1)
  CL__error("Invalid type or size for type_par1 or type_par2");
end
if (type_par1 <> "vinf")
  CL__error("Invalid value for type_par1 (only vinf accepted)");
end

if (typeof(output) <> "string" | size(output, 1) <> 1)
  CL__error("Invalid type or size for output");
end

if (size(output, 2) <> lhs)
  CL__error("Wrong number of output arguments");
end

ParamNames = ["rp", "dinf", "turnang"]; 
if (setdiff(type_par2, ParamNames) <> [])
  CL__error("Invalid value for type_par2");
end

OutputNames = ["vinf", "vp", "dv", ParamNames]; 
if (setdiff(output, OutputNames))
  CL__error("Invalid value for output");
end

// Check / resize arguments
[par1, par2, N] = CL__checkInputs(par1, 1, par2, 1);

// par 1 = vinf
vinf = par1; 
if (find(par1 < 0) <> [])
  CL__error("Invalid value for 1st input parameter")
end


// ------------------
// Computations
// ------------------
// vinf == 0 => return %nan
I = find(vinf <= 0); 
vinf(I) = %nan; 

// STEP 1: par2 => rp
// -----------------
if (type_par2 == "dinf")
  // Check dinf >= 0
  if (find(par2 < 0) <> [])
    CL__error("Invalid value for input parameter (dinf)")
  end
  rp = flybyPar_vinfd2rp(vinf, par2, mu);
  
elseif (type_par2 == "turnang")
  // Check turnang in [0, %pi]
  if (find(par2 < 0 | par2 >  %pi) <> []) 
    CL__error("Invalid value for input parameter (turnang)"); 
  end
  rp = flybyPar_vinfturnang2rp(vinf, par2, mu);

elseif (type_par2 == "rp")
  // Check rp >= 0
  if (find(par2 < 0) <> []) 
    CL__error("Invalid value for input parameter (rp)"); 
  end
  rp = par2;
end

// STEP 2: rp => all
// -----------------
res = struct();
res.vinf = vinf;
res.rp = rp;
[res.dinf, res.turnang, res.vp, res.dv]  = flybyPar_vinfrp2all(vinf, rp, mu);

// set output arguments
for k = 1 : size(output, 2)
  varargout(k) = res(output(k));
end

endfunction