function o = nma_getOrbitParams(r_vec,v_vec,mu)
%find orbit parameters from the velocity and position vectors 

%function o = getOrbitParams(r_vec,v_vec,mu)
%
%Finds orbit parameters from the velosity and position vectors 
%
%INPUT:
%   r_vec : Vector. The position vector of probe. An array of 3 numbers.
%   v_vec : Vector. The velosity vector of probe. An array of 3 numbers.
%   mu     : Scalar. The gravitional parameter. 
%
%OUTPUT:
%
%   a struct that contains these fields
%
%   energy : Scalar. the orbit Energy.
%   h_vec  : vector. angular momentum. 
%   p      : Scalar. orbit semi-latus rectum, or "parameter"
%   a      : Scalar. orbit semi-major axis length in same units as mu.
%   e_vec  : vector. orbit eccentricity vector.
%   c      : scalar. half the distance between 2 foci.
%   rp     : scalar. distance to perigee 
%   rp_vec : vector. a vector pointing from foci to perigree of same
%            length as rp.
%   ra     : scalar. distance to apegee. 
%   ra_vec : vector. similar to rp_vec, but for the apirgee
%   period : scalar. period of probe.
%
%Example Call:
%
% This is an example how to call this function
%
%   mu=1;
%   v_vec = [1 2 -4.4];   % v0_vec = I + 2J  + 4.4K
%   r_vec = [2 3 4];      % r0_vec = 2I + 3J + 4 K
%
%   orbitProp = nma_getOrbitParams(r_vec,v_vec,mu);
%
% Author: Nasser Abbasi
% Version: 1.0

%change history:
%
% name   date    changes
%-----  ------   ---------
% nma   4/30/03  Initial v1.0
% nma   5/3/02   changed it to return a struct

r      = norm(r_vec);
v      = norm(v_vec);

o.energy  = v^2/2 - mu/r;
o.h_vec   = cross(r_vec,v_vec);
o.p       = norm(o.h_vec)^2/mu;
o.e_vec   = cross(v_vec,o.h_vec)/mu - ( r_vec/norm(r_vec) );
e=norm(o.e_vec);

    if(abs(e-1)<eps)  % parabola
        o.c=Inf;
        o.a=Inf;
    else
        o.a       = - mu/(2*o.energy);        
        o.c       = e*o.a;
    end
  

o.rp      = o.p/(1+e);
o.rp_vec  = o.rp* (o.e_vec/norm(o.e_vec));

if(e<1)
   o.ra      = o.a + (o.a*e);
   o.ra_vec  = -(o.ra)*(o.e_vec/norm(o.e_vec));
   o.period  = 2*pi*sqrt( o.a^3/mu );
else
   o.ra      = Inf;
   o.ra_vec  = [Inf Inf Inf];
   o.period  = Inf;
end