# N-Body inputs from Keplerian Orbits

• A
• Shivi20
Shivi20
Hi,
I've trying to figure out how to convert Keplerian orbits into cartesian orbits to use them as input for my code, however I have unsuccessful in doing so,
I wrote a python script for the orbit transformation, but when I convert the Cartesian values back to Keplerian orbits some of the values are not correct.
Can someone help me as to where I might be going wrong?

Keplerian to Cartesian Code:
Python:
import numpy as np
# Constants
G = 4 * np.pi**2  # Gravitational constant in AU^3/yr^2/M_sun
# Given parameters
e1 = 0.05
e2 = 0.87
a1 = 8.0  # AU
a2 = 510.0  # AU
M1 = 20.0  # Msun
M2 = 25.0  # Msun
M3 = 1000.0  # Msun
chi1 = 1.0  # assumed to be phase angle or anomaly
chi2 = 1.0  # assumed to be phase angle or anomaly
O1 = 0.0  # Longitude of ascending node for binary
O2 = np.pi / 2  # Longitude of ascending node for the third body
# Compute the reduced mass and total mass of the binary
mu1 = G * (M1 + M2)  # Gravitational parameter for binary
mu2 = G * (M1 + M2 + M3)  # Gravitational parameter for the third body around binary
# Function to convert Keplerian elements to Cartesian coordinates
def keplerian_to_cartesian(a, e, i, Omega, omega, chi, mu):
# Solve Kepler's equation for eccentric anomaly
E = chi
for _ in range(10):  # Iterative solution
E = chi + e * np.sin(E)

# True anomaly
sin_nu = np.sqrt(1 - e**2) * np.sin(E) / (1 - e * np.cos(E))
cos_nu = (np.cos(E) - e) / (1 - e * np.cos(E))
nu = np.arctan2(sin_nu, cos_nu)

# Distance
r = a * (1 - e**2) / (1 + e * np.cos(nu))

# Position and velocity in orbital plane
r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0])
v_pqw = np.array([-np.sqrt(mu / a) * np.sin(E), np.sqrt(mu / a) * (e + np.cos(E)), 0])

# Rotation matrices
R3_W = np.array([[np.cos(-Omega), -np.sin(-Omega), 0],
[np.sin(-Omega),  np.cos(-Omega), 0],
[0, 0, 1]])

R1_i = np.array([[1, 0, 0],
[0, np.cos(-i), -np.sin(-i)],
[0, np.sin(-i),  np.cos(-i)]])

R3_w = np.array([[np.cos(-omega), -np.sin(-omega), 0],
[np.sin(-omega),  np.cos(-omega), 0],
[0, 0, 1]])

# Rotation matrix from PQW to ECI
R = np.dot(R3_W, np.dot(R1_i, R3_w))

# Convert position and velocity vectors
r_eci = np.dot(R, r_pqw)
v_eci = np.dot(R, v_pqw)

return r_eci, v_eci
# Compute positions and velocities of the binary components
r1_eci, v1_eci = keplerian_to_cartesian(a1, e1, i1, O1, 0, chi1, mu1)
r2_eci, v2_eci = keplerian_to_cartesian(a1, e1, i1, O1, np.pi, chi1, mu1)
# Compute position and velocity of the third body
r3_eci, v3_eci = keplerian_to_cartesian(a2, e2, 0, O2, 0, chi2, mu2)
# Adjust for the center of mass of the binary
r_cm = (M1 * r1_eci + M2 * r2_eci) / (M1 + M2)
v_cm = (M1 * v1_eci + M2 * v2_eci) / (M1 + M2)
# Convert third body relative to the binary center of mass
r3_eci += r_cm
v3_eci += v_cm
print("Body 1 - Position vector (ECI):", r1_eci)
print("Body 1 - Velocity vector (ECI):", v1_eci)
print("Body 2 - Position vector (ECI):", r2_eci)
print("Body 2 - Velocity vector (ECI):", v2_eci)
print("Body 3 - Position vector (ECI):", r3_eci)
print("Body 3 - Velocity vector (ECI):", v3_eci)
Output:
Body 1 -
Position vector (ECI): [ 3.62765682 -0.51761685 -6.8840837 ]
Body 1 - Velocity vector (ECI): [-12.87552775 -0.618391 -8.2243369 ]
Body 2 - Position vector (ECI): [-3.62765682 0.51761685 6.8840837 ]
Body 2 - Velocity vector (ECI): [12.87552775 0.618391 8.2243369 ]
Body 3 - Position vector (ECI): [242.06896743 578.86597711 0.76489819]
Body 3 - Velocity vector (ECI): [6.87271722 8.74136245 0.91381521]

Python:
Cartesian to Keplerian
import numpy as np
def cartesian_to_keplerian(r, v, mu):
"""
Convert Cartesian coordinates to Keplerian elements.

Parameters:
r   : Position vector in Cartesian coordinates (numpy array)
v   : Velocity vector in Cartesian coordinates (numpy array)
mu  : Gravitational parameter (mu = G * M, in AU^3/yr^2)

Returns:
a       : Semi-major axis (in AU)
e       : Eccentricity (dimensionless)
Omega   : Longitude of the ascending node (in radians)
omega   : Argument of periapsis (in radians)
M       : Mean anomaly (in radians)
"""
# Compute the specific angular momentum vector and magnitude
h = np.cross(r, v)
h_mag = np.linalg.norm(h)

# Compute the eccentricity vector and magnitude
e_vec = (np.cross(v, h) / mu) - (r / np.linalg.norm(r))
e = np.linalg.norm(e_vec)

# Compute the semi-major axis (a)
r_mag = np.linalg.norm(r)
v_mag = np.linalg.norm(v)
energy = (v_mag**2 / 2) - (mu / r_mag)
a = -mu / (2 * energy)

# Compute the inclination (i)
i = np.arccos(h[2] / h_mag)

# Compute the node line vector (N) and its magnitude
N = np.array([-h[1], h[0], 0])
N_mag = np.linalg.norm(N)

# Compute the longitude of the ascending node (Omega)
if N_mag != 0:
Omega = np.arccos(N[0] / N_mag)
if N[1] < 0:
Omega = 2 * np.pi - Omega
else:
Omega = 0  # Undefined when N_mag is zero

# Compute the argument of periapsis (omega)
if e != 0 and N_mag != 0:
omega = np.arccos(np.dot(N, e_vec) / (N_mag * e))
if e_vec[2] < 0:
omega = 2 * np.pi - omega
else:
omega = 0  # Undefined when e or N_mag is zero

# Compute the true anomaly (nu)
if e != 0:
nu = np.arccos(np.dot(e_vec, r) / (e * r_mag))
if np.dot(r, v) < 0:
nu = 2 * np.pi - nu
else:
nu = 0  # Circular orbit

# Compute the eccentric anomaly (E)
cos_E = (e + np.cos(nu)) / (1 + e * np.cos(nu))
sin_E = np.sqrt(1 - e**2) * np.sin(nu) / (1 + e * np.cos(nu))
E = np.arctan2(sin_E, cos_E)

# Compute the mean anomaly (M)
M = E - e * np.sin(E)

return a, e, i, Omega, omega, M
# Constants
G = 4 * np.pi**2  # Gravitational constant in AU^3/yr^2/M_sun
# Given parameters for bodies
M1 = 20.0  # Msun
M2 = 25.0  # Msun
M3 = 1000.0  # Msun
# Mu values
mu1 = G * (M1 + M2)  # Gravitational parameter for binary (Body 1 and 2)
mu2 = G * (M1 + M2 + M3)  # Gravitational parameter for the third body (Body 3)
# Cartesian coordinates for Body 1
r1_eci = np.array([3.62765682, -0.51761685, -6.8840837])
v1_eci = np.array([-12.87552775, -0.618391, -8.2243369])
# Cartesian coordinates for Body 2
r2_eci = np.array([-3.62765682, 0.51761685, 6.8840837])
v2_eci = np.array([12.87552775, 0.618391, 8.2243369])
# Cartesian coordinates for Body 3
r3_eci = np.array([242.06896743, 578.86597711, 0.76489819])
v3_eci = np.array([6.87271722, 8.74136245, 0.91381521])
# Convert Cartesian coordinates to Keplerian elements for each body
a1, e1, i1, Omega1, omega1, M1 = cartesian_to_keplerian(r1_eci, v1_eci, mu1)
a2, e2, i2, Omega2, omega2, M2 = cartesian_to_keplerian(r2_eci, v2_eci, mu1)  # Same mu for binary
a3, e3, i3, Omega3, omega3, M3 = cartesian_to_keplerian(r3_eci, v3_eci, mu2)
print("Body 1 - Semi-major axis (a):", a1, "AU")
print("Body 1 - Eccentricity (e):", e1)
print("Body 1 - Inclination (i):", np.rad2deg(i1), "degrees")
print("Body 1 - Longitude of the ascending node (Omega):", np.rad2deg(Omega1), "degrees")
print("Body 1 - Argument of periapsis (omega):", np.rad2deg(omega1), "degrees")
print("Body 1 - Mean anomaly (M):", np.rad2deg(M1), "degrees")
print("\nBody 2 - Semi-major axis (a):", a2, "AU")
print("Body 2 - Eccentricity (e):", e2)
print("Body 2 - Inclination (i):", np.rad2deg(i2), "degrees")
print("Body 2 - Longitude of the ascending node (Omega):", np.rad2deg(Omega2), "degrees")
print("Body 2 - Argument of periapsis (omega):", np.rad2deg(omega2), "degrees")
print("Body 2 - Mean anomaly (M):", np.rad2deg(M2), "degrees")
print("\nBody 3 - Semi-major axis (a):", a3, "AU")
print("Body 3 - Eccentricity (e):", e3)
print("Body 3 - Inclination (i):", np.rad2deg(i3), "degrees")
print("Body 3 - Longitude of the ascending node (Omega):", np.rad2deg(Omega3), "degrees")
print("Body 3 - Argument of periapsis (omega):", np.rad2deg(omega3), "degrees")
print("Body 3 - Mean anomaly (M):", np.rad2deg(M3), "degrees")
Output:
Body 1 -
Semi-major axis (a): 8.009610957285123 AU
Body 1 - Eccentricity (e): 0.0897073385555107
Body 1 - Inclination (i): 94.29999998197292 degrees
Body 1 - Longitude of the ascending node (Omega): 180.0 degrees
Body 1 - Argument of periapsis (omega): 164.36973375231568 degrees
Body 1 - Mean anomaly (M): 68.01013414937296 degrees
Body 2 - Semi-major axis (a): 8.009610957285123 AU
Body 2 - Eccentricity (e): 0.0897073385555107
Body 2 - Inclination (i): 94.29999998197292 degrees
Body 2 - Longitude of the ascending node (Omega): 180.0 degrees
Body 2 - Argument of periapsis (omega): 344.36973375231565 degrees
Body 2 - Mean anomaly (M): 68.01013414937296 degrees
Body 3 - Semi-major axis (a): 5875.725418674007 AU
Body 3 - Eccentricity (e): 0.9921561200716326
Body 3 - Inclination (i): 163.11851976707266 degrees
Body 3 - Longitude of the ascending node (Omega): 67.53655968889326 degrees
Body 3 - Argument of periapsis (omega): 210.87719644086638 degrees
Body 3 - Mean anomaly (M): 1.0592958090738662 degrees

[Code tags added by the Mentors]

Last edited by a moderator:
"Debug my code for me!" is a big ask.
"Oh, and I won't put it in code tags" is a bigger ask.
"And I won't point you to why I think its problematic - just dump the output, You figure it out." is a bigger ask still.
"And its too much bother for me to explain the ideas behind my code. With enough effort, you guys can figure it out from the comments." is again, you guessed it, a big ask.

Welcome to PF.

Shivi20 said:
Can someone help me as to where I might be going wrong?
Can you give us examples of how this code is not working for you?

I'm afraid the code tags include things that aren't code - the outputs.

berkeman
I'm afraid the code tags include things that aren't code - the outputs.
Oops, sorry about that. Working on fixing that now...

berkeman said:
Working on fixing that now...
Why?

I mean, it's nice for you to do it, but it's not your job to fix it. It was the OP's job to write an answerable question, not a tangled mess. I admit I am in a bad mood at the moment, but it shouldn't be the job of PF to fix a bad question.

berkeman
Now that @berkeman has spent the effort to clean this up, I can try and understand it. First, when you write numbers down like 8.00961095728512 AU, that's innumerate nonsense. You are telling us you know where a planet is to within 1/2000 of an inch? Really?

Your variable names make it impossible to debug: M3 is the mass of object 3, but the eccentricity of that object is e2?

My strong advice is to start over, this time with sensible variable names. That alone might tell you where the errors are.

Addendum: you have two variables named e and E. Does thios strike you as a good idea? Especially as E is also really chi.

Last edited:
berkeman

Replies
3
Views
4K
Replies
1
Views
4K
Replies
4
Views
3K
Replies
7
Views
1K
Replies
2
Views
4K
Replies
6
Views
7K
Replies
1
Views
2K
Replies
2
Views
3K
Replies
1
Views
3K
Replies
7
Views
8K