//  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_dvHohmannG(ai,ei,af,ef,posman1,rotation,mu,res)
// Generalized Hohmann transfer
//
// Calling Sequence
// [deltav,dv1,dv2,anv1,anv2] = CL_man_dvHohmannG(ai,ei,af,ef [,posman1,rotation,mu])
// man = CL_man_dvHohmannG(ai,ei,af,ef [,posman1,rotation,mu], res="s")
//
// Description
// <itemizedlist><listitem>
// <p>Computes the maneuvers of a generalized Hohmann transfer. </p>
// <p>This function computes the 2 velocity increments that transform an initial orbit defined 
// by semi-major axis and eccentricity into a final orbit also defined by semi-major axis 
// and eccentricity.</p> 
// <p>The first velocity increment can be at the periapsis or apoapsis of the initial orbit, depending 
// on the value of <b>posman1</b>. </p>
// <p>The second velocity increment is 180 degrees apart.</p>
// <p><b>deltav</b> is the sum of the norms of the velocity increments. </p>
// <p>The velocity increments are given in the "qsw" local orbital frame. </p>
// <p><b>anv1</b> is the true anomaly (on the initial orbit) at the location of the 1st velocity increment, 
// and <b>anv2</b> is the anomaly (on the intermediate orbit) at the location of the 2nd velocity increment.</p>
// <p><b>posman1</b> is used to define the location of the first maneuver. </p>
// <p><b>rotation</b> defines whether the final orbit should have its line of apsides rotated 180deg relative to the initial orbit. </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>
// <para/><inlinemediaobject><imageobject><imagedata fileref="hohmanng.gif"/></imageobject></inlinemediaobject></listitem>
// </itemizedlist>
//
// Parameters
// ai : Semi-major axis of initial orbit [m] (1xN or 1x1)
// ei : Eccentricity of initial orbit (1xN or 1x1)
// af : Semi-major axis of final orbit [m] (1xN or 1x1)
// ef : Eccentricity of final orbit (1xN or 1x1)
// posman1: (optional) Flag indicating the location of the first maneuver: 0 or "per" -> periapsis, 1 or "apo" -> apoapsis. Default is at the periapsis. (1xN or 1x1)
// rotation: (optional) Flag indicating whether the final orbit should have its perigee rotated 180deg relative to the initial orbit (0->no rotation, 1->rotation ; default is 0) (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 : Sum of norms of velocity increments = |dv1|+|dv2| [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)
// anv1 : True anomaly at the location of the first velocity increment (in the initial orbit). [rad] (1xN)
// anv2 : True anomaly at the location of the second velocity increment (in the intermediate orbit). [rad] (1xN)
// man : Structure containing all the output data. 
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_man_dvHohmann
// CL_man_dvBiElliptic
// CL_man_sma
//
// Examples
// // Maneuver at apogee, no rotation
// ai = 7200.e3;
// af = 7000.e3;
// ei = 0.1;
// ef = 0.1;
// [deltav,dv1,dv2,anv1,anv2] = CL_man_dvHohmannG(ai,ei,af,ef,posman1="apo",rotation=0)
//
// // Check results:
// kep = [ai ; ei ; %pi/2 ; 0 ; 0 ; anv1];
// kep1 = CL_man_applyDvKep(kep,dv1);
// kep1(6) = anv2;
// kep2 = CL_man_applyDvKep(kep1,dv2)
//
// // Maneuver at apogee, with rotation
// ai = 7200.e3;
// af = 7000.e3;
// ei = 0.1;
// ef = 0.1;
// [deltav,dv1,dv2,anv1,anv2] = CL_man_dvHohmannG(ai,ei,af,ef,posman1="apo",rotation=1)
//
// // Check results:
// kep = [ai ; ei ; %pi/2 ; 0 ; 0 ; anv1];
// kep1 = CL_man_applyDvKep(kep,dv1);
// kep1(6) = anv2;
// kep2 = CL_man_applyDvKep(kep1,dv2)


// Declarations:

// Code:
if (~exists("posman1", "local")); posman1 = 0; end
if (~exists("rotation", "local")); rotation = 0; end
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

// convert posman type to "real"
if (typeof(posman1) == "string")
  str = posman1; 
  posman1 = %nan * ones(str); 
  posman1(find(str == "per")) = 0; 
  posman1(find(str == "apo")) = 1; 
end

// checks arguments sizes are OK / resizes
[ai,ei,af,ef,posman1,rotation] = CL__checkInputs(ai,1,ei,1,af,1,ef,1,posman1,1,rotation,1); 

if (find(ai <= 0 | ei < 0 | ei >= 1 | af <= 0 | ef < 0 | ef >= 1) <> [])
  CL__error("Invalid input arguments (orbital elements)"); 
end

if (find(posman1 <> 0 & posman1 <> 1) <> [])
  CL__error("Invalid value for ''posman1''"); 
end

if (find(rotation <> 0 & rotation <> 1) <> [])
  CL__error("Invalid value for ''rotation''"); 
end


// sgn1: (1st man at periapsis) => 1, apoapsis => -1
sgn1 = ones(posman1); 
I = find(posman1 == 1); 
sgn1(I) = -1;

// sgnrot: (rotation == 0) => 1, (rotation == 1) => -1 
sgnrot = ones(rotation); 
I = find(rotation == 1); 
sgnrot(I) = -1; 

// 1st man : 
// initial orbit: r1 / r1_opp (maneuver at r1)
// maneuver: r1_opp -> r2_opp (opposite man position)
//
// 2nd man : 
// initial orbit: r2_opp / r1 (maneuver at r2_opp)
// maneuver: r1 -> r2 (opposite man position)

r1     = ai .* (1 - sgn1 .* ei); 
r1_opp = ai .* (1 + sgn1 .* ei); 
r2_opp = af .* (1 + sgn1 .* sgnrot .* ef); 
r2     = af .* (1 - sgn1 .* sgnrot .* ef); 

[dv1, anv1] = CL__man_dvRaps(r1, r1_opp, r2_opp, mu); 
[dv2, anv2] = CL__man_dvRaps(r2_opp, r1, r2, mu); 

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

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


endfunction