//  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 [sma,ecc] = CL_op_rpvp2ae(rp,vp, mu)
// Periapsis radius and velocity at periapsis to semi major axis and eccentricity
//
// Calling Sequence
// [sma,ecc] = CL_op_rpvp2ae(rp,vp [,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes semi major axis and eccentricity from
// periapsis radius and velocity at periapsis (elliptic and hyperbolic orbits).</p>
// <p></p></listitem>
// <listitem>
// <p><b>Note:</b></p>
// <p>The results are correct even if (rp,vp) are data that actually correspond to the apoapsis.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// rp: Periapsis radius [m] (1xN or 1x1)
// vp: Velocity at periapsis [m/s] (1xN or 1x1)
// mu : (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// sma: Semi major axis [m] (1xN)
// ecc: Eccentricity (1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_op_rava2ae
// CL_op_rarp2ae
// CL_op_rpvinf2ae
//
// Examples
// // Elliptic orbit
// rp = 7078.e3;
// vp = 8.e3;
// [sma,ecc] = CL_op_rpvp2ae(rp,vp)
//
// // Hyperbolic orbit
// rp = 7078.e3;
// vp = 12.e3;
// [sma,ecc] = CL_op_rpvp2ae(rp,vp)

// Declarations:

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

Nr = size(rp,1); 
[rp, vp, N] = CL__checkInputs(rp, Nr, vp, Nr); 

if (find(rp <= 0 | vp <= 0) <> []); CL__error("Positive arguments expected"); end;

sma = zeros(1,N); 
ecc = zeros(1,N); 

energy = vp.^2/2 - mu ./ rp;  // calculation of energy to decide whether it is hyperb or ellipse

if (find(energy == 0) <> [])
  CL__error('Invalid input values (parabola case not handled)'); 
end

I = find(energy < 0);   // ellipse
sma(I) = (mu .* rp(I)) ./ (2*mu - rp(I) .* vp(I).^2); 
ecc(I) = abs(sma .* vp(I).^2 - mu) ./ (mu + sma(I) .* vp(I).^2); 

I = find(energy > 0);  // hyperbola
sma(I) = (mu .* rp(I)) ./ (rp(I) .* vp(I).^2 - 2*mu); 
ecc(I) = (rp(I) .* vp(I).^2 - mu) ./ mu; 

endfunction