I can’t get this rocket into orbit. Could somebody help. Please don’t change the initial conditions but the rest is free game. It needs to do this: Update your Finite Difference Method code such that
b) The simulation takes into account atmospheric drag.
c) The spaceship thrusts in the same direction as the velocity.
d) The spaceship has at least two stages, with different specific impulses between 2.5 km/s and 4.5 km/s.
e) The spaceship reaches a stable orbit without lithobraking (make sure to simulate a full orbit to verify!), and stops accelerating.
Here is the code:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
alt = 2 # altitude km
v = 1 # norm of velocity km/s
dt = 0.000001 # theta
Cd = 4 # Drag coefficient
A = 0.0000012 # area
rho = 1.2e+9 # Atmosphere density in kg/km^3
initial_velocity = np.array([v, 0.0, -dt * v])
initial_position = np.array([6378 + alt, 0.0, 0.0]) # (x, y, z) in km
t_span = (0, 1e5) # Time span
gravitational_constant = 398601
mass = 1570000
def dynamics(t, state):
x, y, z, vx, vy, vz, mass = state
r = np.linalg.norm([x, y, z])
gravity_acc = -gravitational_constant * np.array([x, y, z]) / r ** 3
V = np.linalg.norm([vx, vy, vz])
rho = compute_density(z, y) # Updated to compute density based on altitude
if np.allclose([x, y, z], [6378 + alt, 100000.0, 0.0]):
initial_position = np.array([6378 + alt, 0.0, 0.0])
if V != 0:
drag_acc = -0.5 * Cd * A * rho * V ** 2 * np.array([vx, vy, vz]) / V
else:
drag_acc = np.zeros(3)
isp = specific_impulse(V)
rho = compute_density(z, y)
thrust_acc = thrust_acceleration(mass, vx, vy, vz, isp, t)
acceleration = gravity_acc + (drag_acc + thrust_acc) / mass
mdot = np.linalg.norm(thrust_acc) / mass
return [vx, vy, vz, acceleration[0], acceleration[1], acceleration[2], mdot]
def compute_density(z, y):
alt = np.linalg.norm([6378 + z, y, 0.0) - 6378
rho0 = 1.2e+9 # Atmosphere density at sea level in kg/km^3
scale_height = 8.5 # Atmospheric scale height in km
return rho0 / np.exp(alt / scale_height)
def specific_impulse(V):
if 2.5 <= V <= 3:
return 400 # Specific impulse for stage 1 (km/s)
elif 3 < V <= 3.5:
return 200 # Specific impulse for stage 2 (km/s)
elif 3.5 < V <= 4:
return 100 # Specific impulse for stage 3 (km/s)
elif 4 < V <= 4.5:
return 0 # Specific impulse for stage 4 (km/s)
def thrust_acceleration(mass, vx, vy, vz, isp, t):
direction = np.array([vx, vy, vz])
norm = np.linalg.norm(direction)
if norm < 1e-6:
direction = np.array([0, 0, 100])
else:
direction /= norm
thrust_strength = 0.017 * mass
if t >= 7000:
thrust_strength = 0.0 * mass
return thrust_strength * direction
sol = solve_ivp(dynamics, t_span, np.concatenate((initial_position, initial_velocity, [mass])), method='RK45')
plt.figure()
alt = np.sqrt(sol.y[0]**2+ sol.y[1]**2+ sol.y[2]**2)-6378
plt.plot(sol.t, alt)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2], label='Orbit')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.set_title('3-Dimensional Orbit')
ax.legend()
plt.show
()
John G is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.