I am trying to convert between keplerian orbital elements to cartesian position and velocity but something is going wrong. The issue is that I (as a check) compute some keplerian orbital elements with these derived cartiesian position and velocity but the values do not match (I give an example below where the eccentricities do not match which is referred to as check1 in the code). I am trying to code the equation given in,Keplerian_to_cartesian

def kep_2_cart(a,e,i,omega_AP,omega_LAN, EA):

nu = 2*np.arctan(np.sqrt((1+e)/(1-e)) * np.tan(EA/2))
r = a*(1 - e*np.cos(EA))
h =  np.sqrt(mu*a * (1 - e**2))
Om = omega_LAN
w =  omega_AP

X = r*(np.cos(Om)*np.cos(w+nu) - np.sin(Om)*np.sin(w+nu)*np.cos(i))
Y = r*(np.sin(Om)*np.cos(w+nu) + np.cos(Om)*np.sin(w+nu)*np.cos(i))
Z = r*(np.sin(i)*np.sin(w+nu))

p = a*(1-e**2)

V_X = (X*h*e/(r*p))*np.sin(nu) - (h/r)*(np.cos(Om)*np.sin(w+nu) + \
V_Y = (Y*h*e/(r*p))*np.sin(nu) - (h/r)*(np.sin(Om)*np.sin(w+nu) - \
V_Z = (Z*h*e/(r*p))*np.sin(nu) - (h/r)*(np.cos(w+nu)*np.sin(i))

r_vec = np.array([X,Y,Z])
v_vec = np.array([V_X,V_Y,V_Z])

ang_mom_vec = np.cross(r_vec,v_vec)

e_vec = (np.cross(v_vec,ang_mom_vec) - (mu * r_vec/(np.linalg.norm(r_vec))))/mu

#As a check make sure the magitdue of e_vec = e and ang_mom_vec = h
check1 = np.linalg.norm(e_vec) - e
check2 = np.linalg.norm(ang_mom_vec) - h

return [X,Y,Z],[V_X,V_Y,V_Z]
  • 1
    $\begingroup$ @uhoh I have edited the question where I have added a link to the equation I am trying to program. $\endgroup$ – Warrenmovic Jun 22 '20 at 11:08
  • 1
    $\begingroup$ @uhoh sorry about that,I have added a pdf with the equations $\endgroup$ – Warrenmovic Jun 22 '20 at 11:23
  • 1
    $\begingroup$ @uhoh no problem, I really have tried to write out the steps in Python and everything seems ok (in my head), but when I try to get the orbital elements back from the cartesian vectors (outputted) the magnitudes don't match with what I put in...so I am really confused $\endgroup$ – Warrenmovic Jun 22 '20 at 13:44
  • 1
    $\begingroup$ @uhoh not sure if this helpful, but there supposedly a matlab code that does this, uk.mathworks.com/matlabcentral/fileexchange/… $\endgroup$ – Warrenmovic Jun 22 '20 at 13:51
  • 1
    $\begingroup$ Why are you not using the suggested arctan2 function (np.arctan2 in Python+numpy) to ensure that the result is in the correct quadrant, regardless of inputs ? assume this is Python3+ and you are not getting integer division by default? $\endgroup$ – astrosnapper Jun 22 '20 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.