//  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 [isvalid,orbit_type] = CL__oe_isValid(oe_type,oe,mu)
// Internal function to check the validity of orbital elements, also returns the type of orbit
//
// Calling Sequence
// [isvalid,orbit_type] = CL__oe_isValid(oe_type,oe,mu)
//
// Description
// <itemizedlist><listitem>
// <p>The function checks the validity of the given orbital elements.</p>
// <p></p>
// <p>Orbital elements that can be checked are "car", "kep", "cir", "cireq" and "equin"</p>
// <p></p>
// <p>The checks done for cartesian orbital elements are :</p>
// <p>- the given input must be of correct size: 6xN (NB: [] is considered to be invalid!)</p>
// <p>- pos, vel must not be zero vectors</p>
// <p>- pos, vel must not be colinear</p>
// <p></p>
// <p>The checks done for other orbital elements are :</p>
// <p>- the given input must be of correct size: 6xN (NB: [] is considered to be invalid!)</p>
// <p>- sma must be strictly positive</p>
// <p>- e must be positive</p>
// <p>- inc must be in [0,pi]</p>
// <p></p>
// <p>The function additionaly returns the type of orbit:</p>
// <p>1 for elliptical orbits, 2 for hyperbolic and 3 for parabolic orbits</p>
// <p></p>
// <p>Notes: </p>
// <p>- for orbital elements "car", oe = [pos;vel]</p>
// <p>- mu is only needed for cartesian orbital elements !</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// oe_type: (string) Type of orbital elements: "car", "kep", "cir", "cireq" or "equin" (1x1)
// oe: Orbital elements (6xN)
// mu : Gravitational constant. (Can be omitted if oe <> "car") [m^3/s^2]
// isvalid: (boolean) %t if ALL the orbital elements are valid, %f otherwise (1x1)
// orbit_type: Type of orbit: 1 for elliptical orbits, 2 for hyperbolic and 3 for parabolic orbits (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// kep = [7000.e3; 1e-4; 1; %pi/4; 2; %pi/4]
// [isvalid,orbit_type] = CL__oe_isValid("kep",kep);

// Declarations:
EPS_ORB = CL__dataGetEnv("epsOrb", internal=%t);

// Code
if (oe_type <> "car" & oe_type <> "kep" & oe_type <> "cir" &  ...
    oe_type <> "cireq" & oe_type <> "equin")
  CL__error("Invalid type of orbital elements"); 
end

isvalid = %t;
N = size(oe,2);
orbit_type = ones(1,N);

// Check #1 : size
if (size(oe,1) <> 6 & size(oe,1) <> 0)
  CL__error("Invalid size of input argument (number of rows)"); 
end

if (oe == [])
  isvalid = %f; 
  return; // <= RETURN 
end

// Special computation in case of cartesian orbital elements
if (oe_type == "car")
  // Check validity of position, velocity
  r = CL_norm(oe(1:3,:));
  V = CL_norm(oe(4:6,:));
  C = CL_norm(CL_cross(oe(1:3,:),oe(4:6,:)));

  if (find(r == 0 | V == 0 | C == 0) <> [])
    isvalid = %f; 
  end

  // Energy
  E = V.^2 / 2 - mu ./ r;
  
  // ecc = sqrt( 1 + 2.c^2.E/mu^2 )
  // (real just in case)
  ecc = real(sqrt(1 + 2 * C.^2 .* E / mu^2));
end

// Check #2 : sma
if (oe_type <> "car")
  if (find(oe(1,:) <= 0) <> []); 
    isvalid = %f; 
  end
end

// Check #3 : eccentricity
if (oe_type == "kep")
  ecc = oe(2,:);
elseif (oe_type <> "car")
  ecc = sqrt(oe(2,:).^2 + oe(3,:).^2);
end

if (find(ecc < 0) <> []); 
  isvalid = %f; 
end
  
// Parabolic orbits: ecc strictly equal to 1
orbit_type(find(ecc == 1)) = 3;
  
// Hyperbolic orbits: ecc > 1
orbit_type(find(ecc > 1)) = 2;

// Elliptic orbits: ecc < 1 
orbit_type(find(ecc < 1)) = 1;


// Check #4 : inclination
if (oe_type == "kep")
  if (find(oe(3,:) < 0 | oe(3,:) > %pi) <> []); isvalid = %f; end;
elseif (oe_type == "cir")
  if (find(oe(4,:) < 0 | oe(4,:) > %pi) <> []); isvalid = %f; end;
elseif (oe_type == "cireq")
  s = oe(4,:).^2+oe(5,:).^2; // sin(i/2)^2
  if (find(s > 1) <> []); isvalid = %f; end;
elseif (oe_type == "equin")
  // inclination is always in [0,pi] for equinoctial orbits! 
end

  

endfunction