simulate and animate a rocket launch to orbit in 3d in python

To simulate and animate a rocket launch to orbit in 3d, you can use Python along with various libraries like numpy, matplotlib, and mayavi. Here are the general steps you can follow:

  1. Define the rocket's initial position, velocity, and acceleration, and its initial fuel level and thrust.
  2. Define the gravitational force acting on the rocket, which is proportional to the rocket's mass and inversely proportional to the distance between the rocket and the center of the Earth.
  3. Use numerical integration (e.g. Euler's method or Runge-Kutta methods) to update the rocket's position, velocity, acceleration, and fuel level over time.
  4. Use mayavi to draw a 3d plot of the rocket's trajectory in space, along with the Earth and any other celestial bodies involved.
  5. Use matplotlib.animation to create a dynamic, animated plot of the rocket's motion over time, including labels and annotations.

Here's some sample code to get you started:

main.py
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
from mayavi import mlab

# Constants
G = 6.674e-11 # gravitational constant
M_EARTH = 5.972e24 # mass of Earth
R_EARTH = 6371e3 # radius of Earth

# Initial conditions
x0 = 0 # initial x position
y0 = R_EARTH # initial y position (on Earth's surface)
z0 = 0 # initial z position
vx0 = 0 # initial x velocity
vy0 = 0 # initial y velocity
vz0 = 0 # initial z velocity
m0 = 1000 # initial mass (fuel + rocket)
thrust = 1e6 # initial thrust

# Time step
dt = 0.1 # seconds

# Simulation time
tmax = 2000 # seconds
nsteps = int(tmax / dt) + 1

# Arrays for storing data
t = np.zeros(nsteps)
x = np.zeros(nsteps)
y = np.zeros(nsteps)
z = np.zeros(nsteps)
vx = np.zeros(nsteps)
vy = np.zeros(nsteps)
vz = np.zeros(nsteps)
m = np.zeros(nsteps)

# Set initial conditions
t[0] = 0
x[0] = x0
y[0] = y0
z[0] = z0
vx[0] = vx0
vy[0] = vy0
vz[0] = vz0
m[0] = m0

# Simulation loop
for i in range(nsteps-1):

    # Compute gravitational force
    r = np.sqrt(x[i]**2 + y[i]**2 + z[i]**2)
    F_grav = -G * M_EARTH * m[i] / r**2
    Fx_grav = F_grav * x[i] / r
    Fy_grav = F_grav * y[i] / r
    Fz_grav = F_grav * z[i] / r

    # Compute thrust force
    Fx_thrust = 0
    Fy_thrust = 0
    Fz_thrust = 0
    if m[i] > 0:
        Fx_thrust = thrust * vx[i] / np.sqrt(vx[i]**2 + vy[i]**2 + vz[i]**2)
        Fy_thrust = thrust * vy[i] / np.sqrt(vx[i]**2 + vy[i]**2 + vz[i]**2)
        Fz_thrust = thrust * vz[i] / np.sqrt(vx[i]**2 + vy[i]**2 + vz[i]**2)

    # Compute total force
    Fx = Fx_grav + Fx_thrust
    Fy = Fy_grav + Fy_thrust
    Fz = Fz_grav + Fz_thrust

    # Update velocity
    vx[i+1] = vx[i] + Fx / m[i] * dt
    vy[i+1] = vy[i] + Fy / m[i] * dt
    vz[i+1] = vz[i] + Fz / m[i] * dt

    # Update position
    x[i+1] = x[i] + vx[i+1] * dt
    y[i+1] = y[i] + vy[i+1] * dt
    z[i+1] = z[i] + vz[i+1] * dt

    # Update mass
    m[i+1] = m[i] - thrust * dt / np.sqrt(vx[i]**2 + vy[i]**2 + vz[i]**2)

    # Update time
    t[i+1] = t[i] + dt

# 3d plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x/R_EARTH, y/R_EARTH, z/R_EARTH, 'b')
ax.set_xlabel('x (R_EARTH)')
ax.set_ylabel('y (R_EARTH)')
ax.set_zlabel('z (R_EARTH)')
ax.set_title('Rocket trajectory from Earth to orbit')
plt.show()

# Animation
fig = mlab.figure()
mlab.plot3d(x/R_EARTH, y/R_EARTH, z/R_EARTH, t/tmax, tube_radius=0.025)
mlab.show()
2460 chars
103 lines

gistlibby LogSnag