//  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 [M,N] = CL_op_gauss(type_oe,oe,frame, mu)
// Gauss equations
//
// Calling Sequence
// [M, N] = CL_op_gauss(type_oe,oe,frame [,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes Gauss equations for a given type of orbital elements, in the "qsw" or "tnw" frame.</p>
// <p>The function returns M and N such that d(oe)/dt = M*acc + N, where
//    acc = acceleration with components in "qsw", "tnw" or "inertial" frame.</p>
// <p>The "inertial frame" option is only present for convenience (and is less efficient). </p>
// <p></p>
// <p>Available types of orbital elements are: "kep", "cir", "cireq" and "equin".</p>
// <p>See <link linkend="Orbital elements">Orbital elements</link> for the definition of orbital elements</p>
// <p></p>
// <p>Notes: </p>
// <p>- For orbital elements that have singularities (e.g Keplerian orbital elements with an eccentricity of zero)
// the parameters that cannot be computed are set to %nan.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// type_oe: (string) Type of input orbital elements ("kep", "cir", "cireq", "equin" or "car") (1x1)
// oe: Input orbital elements (6xN)
// frame: (string) Coordinates frame for acceleration : "qsw", "tnw" or "inertial". (1x1)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// M: Effects of acceleration on orbital elements (6x3xN)
// N: Mean motion (6xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Les Equations du Mouvement Orbital Perturbe - Cours de l'Ecole GRGS 2002 - Pierre Exertier, Florent Deleflie
//
// See also
// CL_oe_convert
//
// Examples
// // Keplerian orbital elements 
// kep = [7000.e3;0.01;1.8;0.1;0.2;0.3];
// [M,N] = CL_op_gauss("kep",kep,"qsw");
//
// // d(kep)/dt due to some acceleration:
// acc = [1e-3 ; 1e-4; -2e-3]; // in qsw frame
// dkepdt = M*acc + N;
// 
// // Same orbit, circular orbital elements
// [cir, dcirdkep] = CL_oe_kep2cir(kep);
// [M2,N2] = CL_op_gauss("cir",cir,"qsw");
//
// // Check consistency using jacobian:
// M2 - dcirdkep*M   // --> zero


// Declarations:

// Code:

// Check number of input arguments
if (argn(2) <> 3 & argn(2) <> 4)
  CL__error("Wrong number of input arguments");
end

if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// Check input argument "type_oe"
if (typeof(type_oe) <> "string"); CL__error("Invalid input argument: type_oe"); end;
if (type_oe <> "kep" & type_oe <> "cir" & type_oe <> "cireq" & type_oe <> "equin"); 
  CL__error("Invalid type of orbital elements"); 
end

// Check input argument "oe"
[isvalid,type_orbit] = CL__oe_isValid(type_oe,oe);
if (~isvalid); CL__error("Invalid orbital elements"); end
if (find(type_orbit <> 1) <> []); 
  CL__error("Invalid orbital elements (parabolic or hyperbolic orbit)"); 
end

// Check input argument "frame"
if (typeof(frame) <> "string"); CL__error("Invalid input argument: frame"); end
if (frame <> "qsw" & frame <> "tnw" & frame <> "inertial"); 
  CL__error("Invalid frame"); 
end

// Check input argument "mu"
if (size(mu,"*") <> 1); CL__error("Invalid input argument: mu"); end

// frameg = option for CL__op_gauss. 
// Special case if "inertial" => use "qsw" equations (see below)
frameg = frame;
if (frame == "inertial"); frameg = "qsw"; end

// Computation
if (type_oe == "kep")
  [M,N] = CL__op_gaussKep(oe,frameg, mu);
  
elseif (type_oe == "cir")
  [M,N] = CL__op_gaussCir(oe,frameg, mu);
  
elseif (type_oe == "cireq")
  [M,N] = CL__op_gaussCireq(oe,frameg, mu);
  
elseif (type_oe == "equin")
  [M,N] = CL__op_gaussEquin(oe,frameg, mu);
end

// if frame == inertial : 
// doedt = M * acc_qsw = M * (M_inertial2qsw * acc_inertial)
if (frame == "inertial")
  [posvel] = CL_oe_convert(type_oe,"pv",oe, mu); 
  M = M * CL_fr_qswMat(posvel(1:3,:), posvel(4:6,:)); 
end

endfunction