//  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_dvIncRaanCirc(sma,inci,raani,incf,raanf,posman,mu,res)
// Inclination and RAAN maneuver for circular orbits
//
// Calling Sequence
// [deltav,dv,pso] = CL_man_dvIncRaanCirc(sma,inci,raani,incf,raanf [,posman,mu])
// man = CL_man_dvIncRaanCirc(sma,inci,raani,incf,raanf [,posman,mu], res="s")
//
// Description
// <itemizedlist><listitem>
// <p>Computes the dv needed for an inclination and RAAN change for a circular orbit.</p> 
// <p><b>dv</b> is the velocity increment vector (in cartesian coordinates) in the "qsw" local orbital frame. </p>
// <p><b>pso</b> is the argument of latitude of the maneuver.</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
// sma : Semi-major axis [m] (1xN or 1x1)
// inci : Initial inclination [rad] (1xN or 1x1)
// raani : Initial right ascension of ascending node [rad] (1xN or 1x1)
// incf : Final inclination [rad] (1xN or 1x1)
// raanf : Final right ascension of ascending node [rad] (1xN or 1x1)
// posman : (optional) Flag specifying the position of the maneuver: 1 or "n" -> northern hemisphere, -1 or "s" -> southern hemisphere. Default is in the northern hemisphere. (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)
// pso : Argument of latitude of the maneuver [rad] (1xN)
// man : Structure containing all the output data. 
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_man_dvInc
//
// Examples
// sma = 7200.e3;
// inci = 1;
// raani = 2;
// incf = 1.1;
// raanf = 1.9;
// [deltav,dv,pso] = CL_man_dvIncRaanCirc(sma,inci,raani,incf,raanf)
// kep = [sma;0;inci;0;raani;pso];
// kep1 = CL_man_applyDvKep(kep,dv)

// Declarations:

// Code:
if (~exists("posman", "local")); posman = 1; 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 == "n")) = 1; 
  posman(find(str == "s")) = -1; 
end


// checks arguments sizes are OK / resizes
[sma,inci,raani,incf,raanf] = CL__checkInputs(sma,1,inci,1,raani,1,incf,1,raanf,1); 

if (find(sma <= 0 | inci < 0 | inci > %pi | ...
         incf < 0 | incf > %pi) <> [])
  CL__error("Invalid input arguments"); 
end

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


// pso of orbit intersection
// Note: No need to check the output variable 'intersect' ! 
// (default: pso in [0, %pi])
[pso1, pso2] = CL_gm_intersectPlanes(inci, raani, incf, raanf); 

// check/change hemisphere 
I = find(posman == -1); 
pso1(I) = pso1(I) + %pi * ones(I); 
pso2(I) = pso2(I) + %pi * ones(I); 


// V: norm of velocity
V = sqrt(mu ./ sma);

// "qsw" frame - orbit 1 (X_qsw1 = M1 * X_inertiel)
M1 = CL_rot_angles2matrix([3,1,3], [raani; inci; pso1]); 

// "qsw" frame - orbit 2 (X_qsw2 = M2 * X_inertiel)
M2 = CL_rot_angles2matrix([3,1,3], [raanf; incf; pso2]); 


// difference of velocity vectors 
// coordonnates in "qsw" frame of orbit 1 
// s2 = unit vector "s" - orbit 2, coordinates in "qsw" orbit 1
// Note : component on "q" axis should be 0
// (both velocity vectors are perpendicular to radius vector)
s2 = M1 * (M2' * [0;1;0]); 

dv = [V .* s2(1,:); V .* (s2(2,:)-1); V .* s2(3,:)];
deltav = CL_norm(dv); 

pso = pso1;

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


endfunction