function [a,e,inc,capom,omega,capm,f] = orbelm_xv2orbelm(x,y,z,vx,vy,vz);

%% HEADER: orbelm_xv2orbelm.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%Coordinate System: Ecliptic Mean Equinox;
%                    x-y plane is ecliptic, x axis is line of ascending node
%                    of Earth's orbit at J2000, z-axis is perpendicular to
%                    x-y plane and points in direction of Earth's north pole
%                    at J2000
%%Inputs:
%        x: x component of Rp (units of km or AU)
%        y: y component of Rp (units of km or AU)
%        z: z component of Rp (units of km or AU)
%        vx: vx component of Vp (units of km/s or AU/day)
%        vy: vy component of Vp (units of km/s or AU/day)
%        vz: vz component of Vp (units of km/s or AU/day)
%        NOTE: IF Z is not input, X and VX are assumed to be 
%              (n,3) or (3,n) arrays of the Rp and Vp vectors
%
%%Outputs:
%       a: semi-major axis (AU)
%       e: eccentricty
%       inc: incidence angle [radians]
%       capom: longitude of ascending node [radians]
%       omega: argument of periapse [radians]
%       capm: mean anomaly [radians]
%       f: true anomaly [radians]
%
%%Calls: None
%
%%Required Files: None
%
%%Source:
%      Fitzpatrick, p.70 and Murray & Dermote p.53   
%      Converted from Darin Ragozzine IDL coide,
%      which was converted from Fortran
%
%%Author: Alexander G. Hayes (hayes@gps.caltech.edu)
%
%%Change Log: 
%            09/13/08: Draft Code Escalated to V1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% BEGIN CODE: orbelm_xv2orbelm.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%define output arraya
a = zeros(size(x));
a(:) = nan;
e = a;
inc = a;
capom = a;
omega = a;
capm = a;
f = a;
cape = a;
energy = a;

%define gm for inputs in km and km/sec
G = 6.6725985e-11; %gravitational constant in m^3/kg/s^2
solrad = 6.95508e8; %solar radius in meters
au = 1.4959787066e11; %AU in meters
solmass = 1.9891e30; %solar mass in kg
massunit = (au.^3)./(G.*86400.0.^2);
gm = solmass./massunit;

%define threshold for inclination and eccentricity values
tiny = 1.0e-10;

%reformat the input arrays in required
if ~exist('z','var');
    rin = x;
    vin = y;    
    if size(rin,1)==3;
        x = rin(1,:);
        y = rin(2,:);
        z = rin(3,:);
    elseif size(rin,2)==3;
        x = rin(:,1);
        y = rin(:,2);
        z = rin(:,3);
    else;
        disp(['Incorrect Input Size For r!']);
        return;
    end;
    if size(vin,1)==3;
        vx = vin(1,:);
        vy = vin(2,:);
        vz = vin(3,:);
    elseif size(vin,2)==3;
        vx = vin(:,1);
        vy = vin(:,2);
        vz = vin(:,3);
    else;
        disp(['Incorrect Input Size For v!']);
        return;
    end;
    clear rin vin
    %redefine output arraya
    a = zeros(size(x));
    a(:) = nan;
    e = a;
    inc = a;
    capom = a;
    omega = a;
    capm = a;
    f = a;
    cape = a;
    energy = a;
end;


%%input units are assumed to be in km, but check anyway
test_r = sqrt(x(1).^2 + y(1).^2 + z(1).^2);
if test_r > 100;
    %switch first to meters
    x = x.*1000;
    y = y.*1000;
    z = z.*1000;
    vx = vx.*1000;
    vy = vy.*1000;
    vz = vz.*1000;
    %now switch to units of au and au/day
    x = x./au;
    y = y./au;
    z = z./au;
    vx = vx./(au./86400);
    vy = vy./(au./86400);
    vz = vz./(au./86400);
end;

%compute angular momentum H and thereby the inclination INC
hx = y.*vz - z.*vy;
hy = z.*vx - x.*vz;
hz = x.*vy - y.*vx;
h2 = hx.*hx + hy.*hy + hz.*hz;
h = sqrt(h2);
hz( hz > h) = h( hz > h);
inc = acos(hz./h);

%compute longitude of ascending node CAPOM and the argument of latitude u
fac = sqrt(hx.^2 + hy.^2)./h;
ind = find(fac < tiny);
capom(ind) = 0;
u(ind) = atan2(y(ind),x(ind));
ind2 = find( (inc(ind) - pi) < (10.*tiny) );
u(ind(ind2)) = -u(ind(ind2));
ind = find( fac >= tiny);
capom(ind) = atan2(hx,-hy);
u = atan2(z./sin(inc),x.*cos(capom) + y.*sin(capom));
capom(capom < 0) = capom(capom < 0) + 2.*pi;
u(u < 0) = u(u < 0) + 2.*pi;

%compute the radius R and velocity squared V2 and the dot product RD)TV, the
%energy per unit mas ENERGY
r = sqrt(x.*x+y.*y+z.*z);
v2 = vx.*vx + vy.*vy + vz.*vz;
v = sqrt(v2);
vdotr = x.*vx + y.*vy + z.*vz;
energy = 0.5.*v2 - gm./r;

%determine the type of conic section. Yell if not an ellipse
ind_parab = find( abs(energy.*r./gm) < tiny );
if length(ind_parab) > 0;
    disp(['orbelm_xv2orbelm: Parabolic Orbits Are Preset']);
end;
ind_hyper = find( energy > 0 );
if length(ind_hyper) > 0;
    disp(['orbelm_xv2orbelm: Hyperbolic Orbits Are Preset']);
end;
ind_closed = find( energy < 0);
if (length(ind_parab)+length(ind_hyper)+length(ind_closed)) ~= numel(x);
    disp(['orbelm_xv2orbelm: Unknown Orbit Types Found!']);
end;

%calculate the elliptical orbit elements
if length(ind_closed) > 0;
    a(ind_closed) = -0.5.*gm./energy(ind_closed);
    fac = 1 - h2(ind_closed)./(gm.*a(ind_closed));
    ind = find(fac > tiny);
    e(ind_closed(ind)) = sqrt(fac(ind));
    face = (a(ind_closed(ind)) - r(ind_closed(ind)))./...
        (a(ind_closed(ind)).* e(ind_closed(ind)));
    ind2 = find(face > 1);
    cape(ind_closed(ind(ind2))) = 0;
    ind2 = find( (face > -1) & (fac <= 1) );
    cape(ind_closed(ind(ind2))) = acos( face(ind2) );
    ind2 = find( (face <= -1) );
    cape(ind_closed(ind(ind2))) = pi;
    cape( (vdotr < 0) ) = 2.*pi - cape( (vdotr < 0)) ;
    cw = (cos( cape ) - e)./(1- e.*cos(cape) );
    sw = sqrt( 1 - e.*e).*sin(cape)./(1-e.*cos(cape) );
    w = atan2(sw,cw);
    w(w < 0) = w(w < 0) + 2.*pi;
    ind = find( fac <= tiny);
    e(ind_closed(ind)) = 0;
    w(ind_closed(ind)) = u(ind_closed(ind));
    cape(ind_closed(ind)) = u(ind_closed(ind));
    capm(ind_closed) = cape(ind_closed) - e(ind_closed).*sin( cape(ind_closed) );
    omega(ind_closed) = u(ind_closed) - w(ind_closed);
    ind = find( omega(ind_closed) < 0);
    omega(ind_closed(ind)) = omega(ind_closed(ind)) + 2.*pi;
    omega(ind_closed) = omega(ind_closed) - fix( omega(ind_closed)./(2.*pi) )*2.*pi;
    f(ind_closed) = 2.*atan( sqrt( (1 + e(ind_closed))./(1 - e(ind_closed))) .* ...
        tan(cape(ind_closed)./2) );
end;

%calculate hyperbolic orbital paramters
if length(ind_hyper) > 0;
   a(ind_hyper) = 0.5.*gm./energy(ind_hyper);
   fac = h2./(gm.*a);
   e(ind_hyper) = sqrt(1 + fac(ind_hyper));
   tmpf = (a+r)./(a.*e);
   tmpf(tmpf < 1) = 1;
   capf = log(tmpf + sqrt(tmpf.*tmpf -1 ) );
   capf(vdotr < 0) = -capf(vdotr < 0);
   cw = (e - cosh(capf))./(e.*cosh(capf)-1 );
   cw = sqrt(e.*e -1).*sinh(capf)./(e.*cosh(capf) -1);
   w(ind_hyper) = atan2(sw(ind_hyper),cw(ind_hyper));
   w(w <0) = w(w <0) + 2.*pi;
   ind = find( fac(ind_hyper) <= tiny);
   e(ind_hyper(ind)) = 1;
   tmpf(ind_hyper(ind)) = 0.5.*h2(ind_hyper(ind))./gm;
   w(ind_hyper(ind)) = acos(2.*tmpf(ind_hyper(ind))./r(ind_hyper(ind))-1);
   ind2 = find( vdotr(ind_hyper(ind)) < 0);
   w(ind_hyper(ind(ind2))) = 2.*pi - w(ind_hyper(ind(ind2)));
   tmpf(ind_hyper(ind)) = (a(ind_hyper(ind))+r(ind_hyper(ind)))./...
       (a(ind_hyper(ind)).*e(ind_hyper(ind)));
   capf(ind_hyper(ind)) = log(tmpf(ind_hyper(ind)) + ...
       sqrt(tmpf(ind_hyper(ind)).*tmpf(ind_hyper(ind)) -1 ) );
   capm(ind_hyper) = e(ind_hyper).*sinh(capf(ind_hyper)) - capf(ind_hyper);
   omega(ind_hyper) = u(ind_hyper) - w(ind_hyper);
   omega(omega < 0) = omega(omega < 0) + 2.*pi;
   omega(ind_hyper) = omega(ind_hyper) - fix( omega(ind_hyper)./2./pi).*2.*pi;
end;

%calcualte parabolic orbital parametesr
if length(ind_parab) > 0;
   a(ind_parab) = 0.5.*h2(ind_parab)./gm;
   e(ind_parab) = 1;
   w(ind_parab) = acos(2.*a(ind_parab)./r(ind_parab) -1);
   ind = find( vdotr(ind_parab) < 1);
   w(ind_parab(ind)) = 2.*pi - w(ind_parab(ind));
   tmpf = tan(0.5.*w);
   capm(ind_parab) = tmpf(ind_parab).*(1 + tmpf(ind_parab).*tmpf(ind_parab)./3);
   omega(ind_parab) = u(ind_parab) - w(ind_parab);
   omega(omega < 0) = omega(omega < 0) + 2.*pi;
   omega(ind_parab) = omega(ind_parab) - fix( omega(ind_parab)./(2.*pi)).*2.*pi;
end;


%% END CODE: orbelm_xv2orbelm.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%