N-Body inputs from Keplerian Orbits

Shivi20 · Aug 6, 2024

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
i1 = np.deg2rad(94.3)  # radians
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)
    i       : Inclination (in radians)
    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]

Vanadium 50 · Aug 6, 2024

"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.

berkeman · Aug 6, 2024

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?

Vanadium 50 · Aug 6, 2024

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

berkeman · Aug 6, 2024

Vanadium 50 said:

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

Oops, sorry about that. Working on fixing that now...

Vanadium 50 · Aug 6, 2024

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.

Vanadium 50 · Aug 6, 2024

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.

N-Body inputs from Keplerian Orbits

Similar threads

Hot Threads

Recent Insights