//  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_dvApsidesLine(ai,ei,pomi,pomf,posman,mu,res)
// Delta V required to change the line of apsides
//
// Calling Sequence
// [deltav,dv,anv] = CL_man_dvApsidesLine(ai,ei,pomi,pomf [,posman,mu])
// man = CL_man_dvApsidesLine(ai,ei,pomi,pomf [,posman,mu], res="s")
//
// Description
// <itemizedlist><listitem>
// <p>Computes the velocity increment needed to modify the line of apsides.</p>
// <p> The maneuver consists in making the orbit rotate in the orbit plane. Thus, only the argument of periapsis is changed; the semi major axis or the eccentricity are not. </p>
// <p> The output argument <b>dv</b> is the velocity increment 
// expressed in cartesian coordinates in the "qsw" local frame. </p>
// <p><b>anv</b> is the true anomaly at the position of maneuver.</p>
// <p><b>posman</b> can be used to define where the maneuver occurs. </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>
// </listitem>
// </itemizedlist>
//
// Parameters
// ai : Initial semi major axis [m] (1xN or 1x1)
// ei : Initial eccentricity (1xN or 1x1)
// pomi : Initial argument of periapsis [rad] (1xN or 1x1)
// pomf : Final argument of periapsis [rad] (1xN or 1x1)
// posman : (optional) Flag specifying the maneuver location: 0 or "nper" -> near periapsis, 1 or "napo" -> near apoapsis. Default is close to the periapsis. (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 : Norm of velocity increment. [m/s] (1xN)
// dv : Velocity increment in cartesian coordinates in the "qsw" local frame [m/s] (3xN)
// anv : True anomaly at the position of the maneuver [rad] (1xN)
// man : Structure containing all the output data. 
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_man_dvBiElliptic
// CL_man_dvHohmann
// CL_man_dvHohmannG
// CL_man_dvSma
//
// Examples
// ai = 7200.e3;
// ei = 0.1;
// pomi = 1;
// pomf = 1.1;
// [deltav,dv,anv] = CL_man_dvApsidesLine(ai,ei,pomi,pomf,posman=1)
// // Check results:
// anm = CL_kp_v2M(ei,anv);
// kep = [ai; ei; %pi/2; pomi; 0; anm];
// kep1 = CL_man_applyDvKep(kep,dv)


// Declarations:

// Code:
if (~exists("posman", "local")); posman = 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(posman) == "string")
  str = posman; 
  posman = %nan * ones(str); 
  posman(find(str == "nper")) = 0; 
  posman(find(str == "napo")) = 1; 
end

[ai,ei,pomi,pomf,posman] = CL__checkInputs(ai,1,ei,1,pomi,1,pomf,1,posman,1); 

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

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


// 1) maneuver anomaly

// variation of arg of periapsis (between -pi and pi)
// NB: "modulo" required ! 
dpom = CL_rMod(pomf-pomi, -%pi, %pi);

// maneuver near periapsis (default)
anv = dpom / 2; 

// maneuver near apoapsis 
I = find(posman == 1); 
anv(I) = dpom(I) / 2 + %pi * ones(I); 

// anv [0, 2*pi]
anv = CL_rMod(anv, 0, 2*%pi); 


// 2) velocity increment

// cartesian coordinates of dv in "qsw": [deltav, 0, 0]
// with: deltav = -2 * |v| * sin(slope)
// slope = %pi/2 - angle(velocity vector, radius vector)

vsinslope = ei .* sin(anv) .* sqrt(mu./(ai.*(1-ei.^2))); 

// velocity increment in "qsw" (cartesian coordinates)
dv = [-2 * vsinslope; zeros(anv); zeros(anv)]; 
deltav = CL_norm(dv); 

// output
if (res == "d")
  varargout = list(deltav, dv, anv); 
else
  varargout(1) = struct("deltav",deltav, "dv", dv, "anv", anv); 
end


endfunction