//  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 [dv,ecc,tanoe] = CL_ip_insertionDv(vinf,rph,sma,tanoh,mu)
// Delta-V needed for an insertion around a planet
//
// Calling Sequence
// [dv,ecc,tanoe] = CL_ip_insertionDv(vinf,rph,sma [,tanoh,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the delta-v needed for the transfer from an hyperbolic orbit to an elliptical orbit.</p>
// <p>The insertion maneuver is tangential, i.e. the delta-V is parallel to the velocity vector on the hyperbolic orbit. </p>
// <p>The initial orbit is defined by its velocity at infinity (<b>vinf</b>) and periapsis radius (<b>rph</b>).</p>
// <p>The final orbit is defined by its semi major axis(<b>sma</b>).</p>
// <p>The true anomaly of the maneuver on the initial orbit (<b>tanoh</b>) can optionally be specified. By default tanoh=0 (meaning 'at the periapsis').</p>
// <p>The planet is defined by its gravitational constant mu (default is %CL_mu)</p>
// <p></p>
// <p>Note: If the final orbit is cicular, the periapsis is by convention that of the 
// hyperbolic orbit. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// vinf: Velocity on hyperbolic orbit at infinity [m/s]. (1xN)
// rph: Periapsis radius of hyperbolic orbit [m]. (1xN)
// sma: Semi-major axis of target (elliptical) orbit [m]. (1xN)
// tanoh: (optional) True anomaly of the maneuvre (on the hyperbolic orbit) [rad]. Default value is 0. (1xN)  
// mu: (optional) gravitational constant [m3/s2]. Default value is %CL_mu.
// dv: Norm of the delta-v. (1xN)
// ecc: Eccentricity of the final (elliptical) orbit. (1xN)
// tanoe: True anomaly on the elliptical orbit at the time of the maneuver [rad]. (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Insertion around Earth:
// vinf = [5 , 6] * 1.e3;
// eqRad = 6378.e3;
// rph = [eqRad+250.e3 , eqRad+500.e3];
// sma = [eqRad+350.e3 , eqRad+700.e3];
// [dv, ecc, tanoe] = CL_ip_insertionDv(vinf,rph, sma)
//
// // Insertion around Mars:
// eqRad = CL_dataGet("body.Mars.eqRad");
// mu = CL_dataGet("body.Mars.mu");
//
// vinf = [5 , 6] * 1.e3;
// rph = [eqRad+250.e3 , eqRad+500.e3];
// sma = [eqRad+350.e3 , eqRad+700.e3];
// tanoh = [0, 0.1]; // radians
// [dv,ecc,tanoe] = CL_ip_insertionDv(vinf,rph,sma,tanoh,mu)
//

// Declarations:

// Code:
if (~exists('tanoh','local')); tanoh = 0; end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// check inputs
[vinf,rph,sma,tanoh] = CL__checkInputs(vinf,1,rph,1,sma,1,tanoh,1); 

I = find(sma <= 0);
if (I <> []); 
  CL__error("Invalid sma: must be strictly positive");
end

I = find(vinf <= 0);
if (I <> [])
  CL__error("Invalid vinf: must be strictly positive");
end

I = find(rph <= 0);
if (I <> [])
  CL__error("Invalid rph: must be strictly positive");
end

I = find(tanoh <= -%pi | tanoh >= %pi);
if (I <> [])
  CL__error("Invalid tanoh: must be in ]-pi, pi[");
end
  
[nbarg_out] = argn(1);

// Initial velocity (hyperbolic orbit) :
sma_hyp = mu ./ vinf.^2;
ecc_hyp = 1 + rph ./ sma_hyp;
p_hyp = sma_hyp.*(ecc_hyp.^2-1);
r_hyp = %nan * ones(vinf);

I = find(1 + ecc_hyp .* cos(tanoh) > 0);
if (I <> [])
  r_hyp(I) = p_hyp(I) ./ (1 + ecc_hyp(I) .* cos(tanoh(I)));
end
v_hyp = sqrt(mu * (2 ./ r_hyp + 1 ./ sma_hyp)); // velocity
  
// Final velocity (elliptical orbit) :
r_ell = r_hyp;
I = find(sma > r_ell/2);
v_ell = %nan * ones(vinf); // velocity
if (I <> [])
  v_ell(I) = sqrt(mu * (2 ./ r_ell(I) - 1 ./ sma(I)));
end
dv = v_hyp - v_ell;
  
if (nbarg_out >= 2)
  // We use the fact that the slope of the velocity is the same before/after maneuver. 
  // C_ell ./ (r_ell .* v_ell) = C_hyp ./ (r_hyp .* v_hyp)
  C_hyp = sqrt(mu * p_hyp);
  C_ell = C_hyp .* v_ell ./ v_hyp;
  ecc = real(sqrt( 1 - C_ell.^2 ./ (mu*sma))); // real() in case of numerical errors 
end
  
if (nbarg_out >= 3)
  p_ell = sma .* (1 - ecc.^2);
  sin_gamma = ecc_hyp .* sqrt(mu ./ p_hyp) .* sin(tanoh) ./ v_hyp;

  esintanoe = sin_gamma .* v_ell ./ sqrt(mu ./ p_ell);
  ecostanoe = p_ell ./ r_ell - 1;
  
  tanoe = atan(esintanoe, ecostanoe);
end  
  
endfunction