//  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 res = CL_kp_characteristics(sma,ecc,v, er,mu)
// Various Keplerian characteristics
//
// Calling Sequence
// res = CL_kp_characteristics(sma,ecc,v [,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes many keplerian parameters for elliptic or hyperbolic orbits. Parabolic orbits
// are not handled. </p>
// <p>The computed parameters are : </p>
// <p>- standard quantities: </p>
// <p> semi major axis, eccentricity, apoapsis radius, periapsis radius,
// apoapsis altitude, periapsis altitude, velocity at apoapsis, velocity at periapsis, energy,
// law of areas' constant, period, semi minor axis, conic's p parameter, conic's c parameter (angular
// momentum), mean motion, velocity at infinity on an hyperbola, distance from asymptote to
// focus for an hyperbola, escape velocity at apoapsis, escape velocity at periapsis,
// escape delta-V at apoapsis, escape delta-V at periapsis.</p>
// <p>- position-related quantities : </p>
// <p> true anomaly, mean anomaly, eccentric anomaly, time
// from periapsis, radius, x cartesian coordinate, y cartesian coordinate, slope of velocity vector, velocity.</p>
// <p>The output argument <b>res</b> is a structure that contains all the results. </p>
// <p>Parameters that are irrelevant are given the value %nan. Example: res.period is %nan for hyperbolic orbits. </p>
// <p> <b>res</b> is defined as follows:</p>
// <p>- res.orb_type: type of orbit (1->elliptical, 2->hyperbolic)</p>
// <p>- res.sma: semi major axis [m]</p>
// <p>- res.ecc: eccentricity</p>
// <p>- res.ra: apoapsis radius [m]</p>
// <p>- res.rp: periapsis radius [m]</p>
// <p>- res.ha: apoapsis altitude (= res.ra minus planet (equatorial) radius) [m]</p>
// <p>- res.hp: periapsis altitude (= res.rp minus planet (equatorial) radius) [m]</p>
// <p>- res.va: velocity at apoapsis [m/s]</p>
// <p>- res.vp: velocity at periapsis [m/s]</p>
// <p>- res.vinf: velocity at infinity for an hyperbolic orbit (%nan otherwise) [m/s]</p>
// <p>- res.period: orbit's period [s]</p>
// <p>- res.mm: orbit's mean motion [rad/s]</p>
// <p>- res.energy: orbit's energy [J]</p>
// <p>- res.c_area: C constant for the law of areas [m^2/s]</p>
// <p>- res.l_con: conic's l parameter (semi latus rectum)</p>
// <p>- res.c_con: conic's c parameter (linear eccentricity) [m]</p>
// <p>- res.smb: semi minor axis [m]</p>
// <p>- res.d_asymt: distance asymptote-focus for an hyperbolic orbit (%nan otherwise) [m]</p>
// <p>- res.vlp: escape velocity at periapsis [m/s]</p>
// <p>- res.vla: escape velocity at apoapsis [m/s]</p>
// <p>- res.dvp: escape delta-V at periapsis [m/s]</p>
// <p>- res.dva: escape delta-V at apoapsis [m/s]</p>
// <p>- res.v: true anomaly [rad]</p>
// <p>- res.M: mean anomaly [rad]</p>
// <p>- res.E: eccentric anomaly [rad]</p>
// <p>- res.taninf: true anomaly at infinity (only for hyperbolic orbits, %nan otherwise) [rad]</p>
// <p>- res.tperi: time from periapsis [s]</p>
// <p>- res.r: norm of radius vector [m]</p>
// <p>- res.x: x cartesian coordinate [m]</p>
// <p>- res.y: y cartesian coordinate [m]</p>
// <p>- res.slope: slope of velocity vector [rad]</p>
// <p>- res.vel: velocity [m/s]</p>
// <p>NB: the linear eccentricity (res.c_con) is the distance between the center and the focus
// (or one of the two foci). The latus rectum is the chord parallel to the directrix and passing
// through the focus (or one of the two foci). The semi latus rectum (res.l_con) is half
// of the latus rectum.</p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p> 
// <p> A few names have changed but are still recognized. They should be avoided 
// as could disappear in the future. These names are: </p>
// <p> - k : same as energy </p>
// <p> - C_area : same as c_area </p>
// <p> - d_foy : same as d_asymt</p>
// <p> - anvinf : same as taninf</p>
// <p> - date: same as tperi</p>
// <p> - R: same as r </p>
// <p> - v_slope: same as slope </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// sma: semi major axis [m] (1xN or 1x1)
// ecc: eccentricity  (1xN or 1x1)
// v: true anomaly [rad] (1xN or 1x1)
// mu: (optional) gravitational constant [m^3/s^2] (default value is %CL_mu)
// er: (optional) planet equatorial radius [m] (default is %CL_eqRad)
// res: structure containing all the results (1xN)
//
// See also
// CL_kp_params
// CL_kp_E2v
// CL_kp_v2E
// CL_kp_E2M
// CL_kp_M2E
// CL_op_rarp2ae
// CL_op_rava2ae
// CL_op_rpvinf2ae
// CL_op_rpvp2ae
//
// Authors
// CNES - DCT/SB
//
// Examples
// sma = 7078.e3;
// ecc = 0.0001;
// v = %pi/4;
// res = CL_kp_characteristics(sma,ecc,v)

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

// Code:
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

Nsma = size(sma,2);
Necc = size(ecc,2);
Nv = size(v,2);
N = max(Nsma,Necc,Nv);

if ~( (Nsma==1|Nsma==N) & (Necc==1|Necc==N) & (Nv==1|Nv==N) ) 
   CL__error('Wrong sizes for input arguments'); 
end

if (Nsma==1); sma = sma*ones(1,N); end
if (Necc==1); ecc = ecc*ones(1,N); end
if (Nv==1); v = v*ones(1,N); end


// keplerian parameters
  orb_type = zeros(ecc);
  ih = find(ecc>1);
  ie = find(ecc<1);
  orb_type(ih) = 2; //hyperbolic orbit
  orb_type(ie) = 1; //elliptic orbit

  if find(ecc >= 1-EPS_ORB.parab & ecc <= 1+EPS_ORB.parab) ~= [] 
     CL__error('Parabola case not handled');
  end

  // calculation of main and secondary kepler parameters
  [ra,rp,ha,hp,va,vp,nmoy,period,vinf,anvinf,dinf,vlibp,vliba,l_con,c_con,k,caire,smb,..
  anm,anex,temps,r,x,y,pen,vit] = CL__kp_keplerParams(sma,ecc,v,orb_type,er,mu)

  // calculation of velocity increments (apoapsis and periapsis) only for elliptic orbits
  delta_va = %nan * ones(va)
  delta_vp = %nan * ones(vp)
  delta_va(ie) = vliba(ie)-va(ie)
  delta_vp(ie) = vlibp(ie)-vp(ie)

  res = struct( ..
    'orb_type', orb_type, ..
    'sma', sma, ..
    'ecc', ecc, ..
    'ra', ra, ..
    'rp', rp, ..
    'ha', ha, ..
    'hp', hp, ..
    'va', va, ..
    'vp', vp, ..
    'vinf', vinf, ..
    'period', period, ..
    'mm', nmoy, ..
    'k', k, ..  // DEPRECATED
    'C_area', caire, .. // DEPRECATED
    'l_con', l_con, ..
    'c_con', c_con, ..
    'smb', smb, ..
    'd_foy', dinf, .. // DEPRECATED
    'vlp', vlibp, ..
    'vla', vliba, ..
    'dvp', delta_vp, ..
    'dva', delta_va, ..
    'v', v, ..
    'M', anm, ..
    'E', anex, ..
    'anvinf', anvinf, .. // DEPRECATED
    'date', temps, .. // DEPRECATED
    'R', r, .. // DEPRECATED
    'x', x, ..
    'y', y, ..
    'v_slope', pen, .. // DEPRECATED
    'vel', vit, ..
    'energy', k, ..   // NEW NAME 
    'c_area', caire, ..  // NEW NAME 
    'd_asymt', dinf, .. // NEW NAME 
    'taninf', anvinf, .. // NEW NAME 
    'tperi', temps, .. // NEW NAME 
    'r', r, .. // NEW NAME 
    'slope', pen .. // NEW NAME 
  );


endfunction