//  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_man_dvBiElliptic(ai,af,rt,mu,res)
// Bi-elliptic transfer
//
// Calling Sequence
// [deltav,dv1,dv2,dv3,anv1,anv2,anv3] = CL_man_dvBiElliptic(ai,af,rt [,mu])
// man = CL_man_dvBiElliptic(ai,af,rt [,mu], res="s")
//
// Description
// <itemizedlist><listitem>
// <p>Computes the maneuvers of a bi-elliptical transfer from a circular 
// orbit with semi-major axis <b>ai</b> 
// to a circular orbit with semi-major axis <b>af</b>.</p> 
// <p> The apogee radius of the elliptical transfer orbit is 
// <b>rt</b>. </p>
// <p><b>deltav</b> is the sum of the norms 
// of the velocity increments.</p>
// <p>Velocity increments are expressed in cartesian coordinates in the "qsw" local frame. </p>
// <p>If the argument <b>res</b> is present and is equal to "s", all the output data are returned in a structure. </p>
// <p><inlinemediaobject><imageobject><imagedata fileref="bielliptic.gif"/></imageobject></inlinemediaobject></p>
// </listitem>
// </itemizedlist>
//
// Parameters
// ai : Semi-major axis of initial circular orbit  [m] (1xN or 1x1)
// af : Semi-major axis of final circular orbit [m] (1xN or 1x1)
// rt : Radius at the position of the second maneuver [m] (1xN or 1x1)
// mu : (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// res : (string, optional) Type of output: "d" or "s" for . Default is "d".
// deltav : Som of norms of velocity increments (=|dv1|+|dv2|+|dv3]) [m/s] (1xN)
// dv1 : First velocity increment in cartesian coordinates in the "qsw" frame [m/s] (3xN)
// dv2 : Second velocity increment in cartesian coordinates in the "qsw" frame [m/s] (3xN)
// dv3 : Third velocity increment in cartesian coordinates in the "qsw" frame [m/s] (3xN)
// anv1 : True anomaly at the position of the 1st maneuver. The initial orbit is circular so the value is arbitrary and set to 0  (1xN)
// anv2 : True anomaly at the position of the 2nd maneuver (either 0 or %pi) [rad] (1xN)
// anv3 : True anomaly at the position of the 3rd maneuver (either 0 or %pi) [rad] (1xN)
// man : Structure containing all the output data. 
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Orbital Mechanics for engineering students, H D Curtis, Chapter 6
//
// See also
// CL_man_dvHohmann
// CL_man_dvHohmannG
// CL_man_dvSma
//
// Examples
// // 7000 km to 98 000km through a 280 000 transfer orbit:
// ai = 7000.e3;
// af = 98000.e3;
// rt = 280000.e3;
// [deltav,dv1,dv2,dv3,anv1,anv2,anv3] = CL_man_dvBiElliptic(ai,af,rt)
//
// // Check results:
// kep = [ai ; 0 ; %pi/2 ; 0 ; 0 ; anv1];
// kep1 = CL_man_applyDvKep(kep,dv1);
// kep1(6) = anv2;
// kep2 = CL_man_applyDvKep(kep1,dv2);
// kep2(6) = anv3;
// kep3 = CL_man_applyDvKep(kep2,dv3)
//
// // Same example with a Hohmann transfer:
// // more expensive !
// [deltav,dv1,dv2,anv1,anv2] = CL_man_dvHohmann(ai,af) 

// Declarations:

// Code:
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end
if (~exists("res", "local")); res = "d"; end

// check "res" argument
if (res <> "d" & res <> "s"); 
  CL__error("Invalid value for argument ''res''"); 
end
if (argn(1) > 1 & res == "s"); 
  CL__error("Invalid number of output arguments"); 
end

// checks arguments sizes are OK / resizes
[ai,af,rt] = CL__checkInputs(ai,1,af,1,rt,1); 

if (find(ai <= 0 | af <= 0 | rt <= 0) <> [])
  CL__error("Invalid input arguments"); 
end

// 1st man : 
// initial orbit: ai / ai
// maneuver: ai -> rt (opposite man position)
//
// 2nd man : 
// initial orbit: rt / ai (maneuver at rt)
// maneuver: ai -> af (opposite man position)
//
// 3rd man : 
// initial orbit: af / rt (maneuver at af)
// maneuver: rt -> af (opposite man position)

[dv1, anv1] = CL__man_dvRaps(ai, ai, rt, mu); 
[dv2, anv2] = CL__man_dvRaps(rt, ai, af, mu); 
[dv3, anv3] = CL__man_dvRaps(af, rt, af, mu); 

// sum of dv norms
deltav = CL_norm(dv1) + CL_norm(dv2) + CL_norm(dv3); 

// output
if (res == "d")
  varargout = list(deltav, dv1, dv2, dv3, anv1, anv2, anv3); 
else
  varargout(1) = struct("deltav",deltav, "dv1", dv1, "dv2", dv2, "dv3", dv3, "anv1", anv1, "anv2", anv2, "anv3", anv3); 
end


endfunction