;Edward Burin des Roziers
;9/15/03
;converts from cartesian coordinates to keplerian elements

function c2kr, r_vector, v_vector

mu = double(3.986004e5); Km^3/s^2
pi =double(3.14159265)

v = sqrt(dotp(v_vector, v_vector))
r = sqrt(dotp(r_vector, r_vector))

h_vector = crossp(r_vector, v_vector)
h = sqrt(dotp(h_vector, h_vector))

hv = [0,0,1]

n_vector = crossp(hv, h_vector)
n = sqrt(dotp(n_vector, n_vector))

e_vector = ((v^2-mu/r) * r_vector - dotp(r_vector,v_vector) * v_vector) / mu
e = sqrt(dotp(e_vector, e_vector))

energy = v^2/2 - mu/r

if (e ne 1.0) then begin
    a = -1 * mu / (2 * energy)
    p = a * (1 - e * e)
endif else begin
    p = h * h / mu
    a = 99999e15
endelse

i = acos(h_vector[2] / h)

if ((e lt 1e-8) or (abs(i) lt 1e-8)) then begin
    print, 'Orbit cannot be circular or equatorial'
endif else begin

Om = acos(n_vector[0] / n)

if (n_vector[1] lt 0) then begin
      Om = 2 * pi - Om
endif

w = acos(dotp(n_vector, e_vector) / (n * e))

if (e_vector[2] lt 0) then begin
      w = 2 * pi - w
endif

nu = acos(dotp(e_vector, r_vector) / (e * r))

if (dotp(r_vector,v_vector) lt 0) then begin
      nu = 2 * pi - nu
endif

EA = 2.*atan(((1-e)/(1+e))^.5*tan(nu/2.))

;Tp = t - (sqrt(a^3./mu))*(EA - e*sin(EA))

d = sqrt(a^3/mu)

aeiomwnu = [a,e,i,Om,w,nu]
return, aeiomwnu

endelse

end