- #1
UlrickAKobbeltvedt
- 0
- 4
Hi everyone,
I'm working on a Python simulation that aims to replicate the Apollo 11 translunar injection (TLI) in a 3-dimensional geocentric system. My model includes Earth, the Moon, and a Saturn V–derived Apollo 11 rocket. The simulation starts at the moment of TLI, using initial positions and velocities obtained from NASA Horizons (Earth and Moon data from July 16, 1969, 16:16:16 UT) along with estimated state vectors for the rocket at low Earth orbit just before the injection burn.
My simulation details:
The problem:Despite double-checking the input data and calculations for the gravitational forces, thrust, and mass changes, the simulated trajectory does not match the expected return orbit. Instead of achieving a proper return trajectory after the injection, the rocket’s path seems off.
My questions:
I'm working on a Python simulation that aims to replicate the Apollo 11 translunar injection (TLI) in a 3-dimensional geocentric system. My model includes Earth, the Moon, and a Saturn V–derived Apollo 11 rocket. The simulation starts at the moment of TLI, using initial positions and velocities obtained from NASA Horizons (Earth and Moon data from July 16, 1969, 16:16:16 UT) along with estimated state vectors for the rocket at low Earth orbit just before the injection burn.
My simulation details:
- Dynamics: I model the gravitational forces between the bodies using Newton’s law and incorporate a thrust force during the TLI burn.
- Burn specifics: The burn lasts 350 seconds using a J‑2 engine with a full thrust of 1.033×10⁶ N, a specific impulse (I_sp) of 421 s, and a calculated mass flow (≈250 kg/s). I also update the rocket’s mass during the burn.
- Integration: I employ a 4th-order Runge-Kutta (RK4) method for integrating the equations of motion.
- Data: The Earth and Moon initial positions/velocities are read from CSV files (with Horizons data), and the rocket’s starting state is estimated to be in the same reference frame.
The problem:Despite double-checking the input data and calculations for the gravitational forces, thrust, and mass changes, the simulated trajectory does not match the expected return orbit. Instead of achieving a proper return trajectory after the injection, the rocket’s path seems off.
My questions:
- Force Calculations: Are there any common pitfalls in computing the gravitational or thrust forces (e.g., unit conversion issues) that might lead to significant discrepancies?
- Integration Method: Could the RK4 implementation be a source of error in this context? Any suggestions for validating the integration accuracy?
- Thrust & Mass Update: Is the way I'm applying the thrust (with a fixed desired vector) and updating the mass during burn physically sound for this simulation?
- General Approach: Are there any aspects of the simulation setup (initial conditions, reference frame conversions, etc.) that could be problematic when trying to replicate the Apollo 11 TLI?
Python:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import os
import time
import csv
# Dynamisk sti-håndtering
fil_sti = os.path.join(os.path.dirname(__file__), "CSVresults")
def parse_horizons_initial_data(fil_sti):
if not os.path.exists(fil_sti):
raise FileNotFoundError(f"File not found: {fil_sti}")
with open(fil_sti, "r") as file:
lines = file.readlines()
start_idx = next((i for i, line in enumerate(lines) if "$SOE" in line), None)
if start_idx is None:
raise ValueError("Missing $SOE in the file")
numerical_line = lines[start_idx + 1].strip()
parts = numerical_line.split(",")
if len(parts) < 8:
raise ValueError(f"Line has insufficient data: {numerical_line}")
x, y, z = map(float, parts[2:5])
vx, vy, vz = map(float, parts[5:8])
return [x, y, z], [vx, vy, vz]
# Konstant for gravitasjon (km^3/(kg·s^2))
G = 6.67430e-20
def gravitational_force(body1, body2):
r_vec = body2.position - body1.position
r_mag = np.linalg.norm(r_vec)
if r_mag == 0:
return np.zeros(3)
force_mag = G * body1.mass * body2.mass / r_mag**2
return force_mag * (r_vec / r_mag)
class Body:
def __init__(self, name, mass, position, velocity, radius=1.0):
self.name = name
self.mass = mass # kg
self.position = np.array(position, dtype=float) # km
self.velocity = np.array(velocity, dtype=float) # km/s
self.acceleration = np.zeros(3, dtype=float) # km/s²
self.radius = radius # km
# --- TLI-parametere (oppdatert modell) ---
burn_start_time = 0 # TLI starter med en gang
burn_duration = 350 # Burn varer i 350 sekunder
# Full engine thrust for J‑2 in Newtons
engine_thrust_newton = 1.033e6 # 1.033×10⁶ N
I_sp = 421.0 # spesifikt impuls (s)
g0 = 9.81 # m/s²
# Calculate mass_flow to produce full engine thrust:
# thrust = mass_flow * I_sp * g0 --> mass_flow = thrust / (I_sp * g0)
mass_flow = engine_thrust_newton / (I_sp * g0) # approx. 250 kg/s
# Oppdater startmassen slik at drivstoffmengden stemmer overens med burn_duration
rocket_dry_mass = 65200.0 # kg (etter forbruk av drivstoff, S‑IVB tørrmasse)
rocket_start_mass = rocket_dry_mass + mass_flow * burn_duration # beregnet total masse ved start
# Etter dump av S‑IVB-tørrmasse, nyttig last for videre ferd
final_rocket_mass = 45700.0 # kg
# Definer pre-TLI og ønsket post-TLI hastighetsvektorer
preTLI = np.array([0.3, 7.79, -0.2]) # km/s (før TLI)
postTLI = np.array([-1.5, 9.6, 5.5]) # km/s (ønsket etter TLI)
delta_v_desired = postTLI - preTLI # Ønsket Δv-vektor
desired_thrust_dir = delta_v_desired / np.linalg.norm(delta_v_desired)
# Nå representerer desired_thrust_dir den retningen vi ønsker thrusten skal virke
def compute_accelerations(bodies, current_time):
for body in bodies:
total_force = np.zeros(3)
# Gravitasjonskraft (for alle interaksjoner)
for other_body in bodies:
if body != other_body:
total_force += gravitational_force(body, other_body)
# Legg til thrust for raketten under burn
if body.name == "Rocket" and (burn_start_time <= current_time < burn_start_time + burn_duration):
thrust_dir = desired_thrust_dir
# Bruk full J‑2 thrust
thrust_newton = engine_thrust_newton
thrust = thrust_newton / 1000.0 # konverterer til kg·km/s² (1 N = 0.001 kg·km/s²)
total_force += thrust * thrust_dir
body.acceleration = total_force / body.mass
def rk4_step(bodies, dt, current_time):
init_pos = [body.position.copy() for body in bodies]
init_vel = [body.velocity.copy() for body in bodies]
# K1
compute_accelerations(bodies, current_time)
k1_pos = [dt * body.velocity for body in bodies]
k1_vel = [dt * body.acceleration for body in bodies]
# K2
for i, body in enumerate(bodies):
body.position = init_pos[i] + 0.5 * k1_pos[i]
body.velocity = init_vel[i] + 0.5 * k1_vel[i]
compute_accelerations(bodies, current_time + 0.5 * dt)
k2_pos = [dt * body.velocity for body in bodies]
k2_vel = [dt * body.acceleration for body in bodies]
# K3
for i, body in enumerate(bodies):
body.position = init_pos[i] + 0.5 * k2_pos[i]
body.velocity = init_vel[i] + 0.5 * k2_vel[i]
compute_accelerations(bodies, current_time + 0.5 * dt)
k3_pos = [dt * body.velocity for body in bodies]
k3_vel = [dt * body.acceleration for body in bodies]
# K4
for i, body in enumerate(bodies):
body.position = init_pos[i] + k3_pos[i]
body.velocity = init_vel[i] + k3_vel[i]
compute_accelerations(bodies, current_time + dt)
k4_pos = [dt * body.velocity for body in bodies]
k4_vel = [dt * body.acceleration for body in bodies]
for i, body in enumerate(bodies):
body.position = init_pos[i] + (k1_pos[i] + 2*k2_pos[i] + 2*k3_pos[i] + k4_pos[i]) / 6
body.velocity = init_vel[i] + (k1_vel[i] + 2*k2_vel[i] + 2*k3_vel[i] + k4_vel[i]) / 6
# Masseendring under burn: Reduserer kun drivstoffmassen gradvis
for body in bodies:
if body.name == "Rocket" and (burn_start_time <= current_time < burn_start_time + burn_duration):
body.mass = max(rocket_dry_mass, body.mass - mass_flow * dt)
return current_time + dt
# --- Opprett legemene ---
bodies = []
""" fikk dette som gjør om på fixed earth koordinater til inertial earth koordinater
[6573.207043182418, -154.55746834527565 , 7.4267378947397695]
[-6573.207043182418, 154.55746834527565 , -7.4267378947397695]
"""
# Rakett med oppgitt pre-TLI starttilstand
rocket = Body(
name="Rocket",
mass=rocket_start_mass, # beregnet total masse ved start
position=[6573.207043182418, -154.55746834527565 , 7.4267378947397695], # pre-TLI posisjon (km) ganget med 0.93 for en grunn
velocity=[0.3, 7.79, -0.2], # pre-TLI hastighet (km/s)
radius=5
)
bodies.append(rocket)
# Bruk geosentriske data for Jorda og Månen (fra Horizons-data)
planet_data = {
"Earth": {"file": "Earth1969_results.csv", "mass": 5.972e24},
"Moon": {"file": "Luna1969_results.csv", "mass": 7.348e22},
}
for planet, data in planet_data.items():
file_path = os.path.join(fil_sti, data["file"])
position, velocity = parse_horizons_initial_data(file_path)
bodies.append(Body(planet, data["mass"], position, velocity))
# --- Lister for logging ---
time_list = []
grav_force_list = []
rocket_speed_list = []
min_distance_list = [] # Logg avstand mellom Rocket og Moon per tidssteg
mass_list = [] # Logg rakettens masse over tid
min_distance = np.inf
min_distance_time = None
# Tidsparametere for simuleringen
dt = 5 # tidssteg i sekunder
simulation_time = 2*24*60*60 # simulerer ca. 1 dag
current_time = 0
steps = int(simulation_time / dt)
progress_bar = tqdm(total=steps, desc="Simulation Progress", unit="steps")
# Flag for å sikre at dumpen (separasjonen) bare utføres én gang
separation_done = False
# Lagrer banedata (trajectories)
trajectories = {body.name: {"x": [], "y": [], "z": []} for body in bodies}
for body in bodies:
trajectories[body.name]["x"].append(body.position[0])
trajectories[body.name]["y"].append(body.position[1])
trajectories[body.name]["z"].append(body.position[2])
while current_time < simulation_time:
current_time = rk4_step(bodies, dt, current_time)
progress_bar.update(1)
# Finn referanser til Rocket og Moon
rocket_obj = next((b for b in bodies if b.name == "Rocket"), None)
moon_obj = next((b for b in bodies if b.name == "Moon"), None)
# Etter at burn-perioden er over, utfør dumpen av den resterende S‑IVB-tørrmassen
if current_time >= burn_start_time + burn_duration and not separation_done:
if rocket_obj is not None:
rocket_obj.mass = final_rocket_mass # Sett rakettmassen til 45700 kg
separation_done = True
if rocket_obj is not None and moon_obj is not None:
F_rm = gravitational_force(rocket_obj, moon_obj)
F_rm_mag = np.linalg.norm(F_rm) # i kg·km/s²
rocket_speed = np.linalg.norm(rocket_obj.velocity) # km/s
distance = np.linalg.norm(rocket_obj.position - moon_obj.position)
min_distance_list.append(distance)
if distance < min_distance:
min_distance = distance
min_distance_time = current_time
time_list.append(current_time)
grav_force_list.append(F_rm_mag)
rocket_speed_list.append(rocket_speed)
mass_list.append(rocket_obj.mass) # Logg rakettens masse
for body in bodies:
trajectories[body.name]["x"].append(body.position[0])
trajectories[body.name]["y"].append(body.position[1])
trajectories[body.name]["z"].append(body.position[2])
progress_bar.close()
print("Siste gravitasjonskraft mellom Rocket and Moon (kg·km/s²):", grav_force_list[-1])
print("Siste gravitasjonskraft i Newton:", grav_force_list[-1] * 1000)
print("Last velocity (km/s):", rocket_speed_list[-1])
print("minimal distance between Rocket and Moon (km):", min_distance)
print("Time of minimal distance (s):", min_distance_time)
def plot_trajectories(trajectories):
# Plotter med Jorda i sentrum ved å subtrahere Jordas initielle posisjon
earth_initial_x = trajectories["Earth"]["x"][0]
earth_initial_y = trajectories["Earth"]["y"][0]
fig, ax = plt.subplots(figsize=(10, 10))
earth_radius_km = 6378
circle = plt.Circle((0, 0), earth_radius_km, fill=False, color="black", linestyle="--")
ax.add_patch(circle)
for name, coords in trajectories.items():
x_vals = np.array(coords["x"]) - earth_initial_x
y_vals = np.array(coords["y"]) - earth_initial_y
ax.plot(x_vals, y_vals, label=name, linewidth=1)
ax.set_xlabel("X (km), Earth-Centered")
ax.set_ylabel("Y (km), Earth-Centered")
ax.set_title("2D Orbital Trajectories (Earth at (0,0))")
ax.legend()
ax.grid(True, which='both', linestyle='--', linewidth=0.5)
all_x = []
all_y = []
for name, coords in trajectories.items():
all_x.extend(np.array(coords["x"]) - earth_initial_x)
all_y.extend(np.array(coords["y"]) - earth_initial_y)
margin = 1.1
ax.set_xlim(margin * min(all_x), margin * max(all_x))
ax.set_ylim(margin * min(all_y), margin * max(all_y))
ax.set_aspect("equal", adjustable="datalim")
plt.show()
plot_trajectories(trajectories)Attachments
Last edited by a moderator:
