Use the tabs below to toggle between the output plots and the code (script.py).
# Welcome!
# Below is some code that imports useful packages. It's essentially a toolbox of pre-built code you can use. You don't need to worry about this section though.
# --------imported packages----------
import numpy as np
import matplotlib.pyplot as plt
# -----------------------------------
# You may skip this section
# --------science mumbo jumbo----------
radius_Nalorp = 1 # radius of planet in Nalorp meters
def satellite_trajectory_calc(initial_veloc):
dt = 0.01 # time step between each position calculation
steps = 5000 # total number of steps (defines length of the trajectory)
GM = 1 # gravitational acceleration normalized to 1
initial_pos = np.array([0, 1.05])
r = np.zeros((steps, 2)) # array of positions of the rocket
r[0] = initial_pos # initialize the first position with the function input
dotr = initial_veloc # initialize the first velocity with the function input
# numerical integration:
for i in range(steps - 1):
normR = np.sqrt(np.sum(r[i, :] ** 2))
if normR <= radius_Nalorp and i > 0:
r = r[:i + 1, :]
break
if normR >= 6:
r = r[:i + 1, :]
break
if (np.abs(r[i, 0] - r[0, 0]) < 0.01) and (np.abs(r[i, 1] - r[0, 1]) < 0.01) and i > 1:
r = r[:i + 1, :]
break
ddotr = -GM * r[i] / normR**3
dotr = ddotr * dt + dotr
r[i + 1] = dotr * dt + r[i]
return r, normR
def blast_off_start_traj(fuel, angle):
m_to_normalized_m = 1 / 6e6
fuelkg = fuel * 1000 # metric tons to kg
rocket_weight = 2e6 # (kg)
heat_of_combustion = 4.2e14 # (J/kg)
initVelMagnitude = (2 * heat_of_combustion * fuelkg / rocket_weight) ** 0.5 * m_to_normalized_m
yVel = np.sin(angle * np.pi / 180) * initVelMagnitude
xVel = np.cos(angle * np.pi / 180) * initVelMagnitude
return np.array([xVel, yVel])
def plot_trajectory(init_velocity):
traj, finalRadius = satellite_trajectory_calc(init_velocity)
theta = np.linspace(0, 2 * np.pi, 100)
xNalorp = np.cos(theta)
yNalorp = np.sin(theta)
fig, axs = plt.subplots(1)
axs.plot(xNalorp, yNalorp, 'green')
axs.plot(traj[:, 0], traj[:, 1], 'r:')
axs.set_xlim((-3, 3))
axs.set_ylim((-3, 3))
axs.set_aspect(1)
plt.xticks([], [])
plt.yticks([], [])
[plt.gca().spines[s].set_visible(False) for s in plt.gca().spines]
if launch_pad_activation:
axs.text(-2.9, 2.9, 'Task failed succesfully:\ncrash landed in mountains.\nUse "launch initiated" as your answer to the next lesson.', fontsize=14,
verticalalignment='top')
else:
axs.text(-2.9, 2.9,
'Launch protocol aborted:\nlaunch pad not activated.',
fontsize=14,
verticalalignment='top')
plt.gca().clear()
# --------end of science mumbo jumbo--------
# --------Important stuff--------
# To activate launch pad set: False -> True
launch_pad_activation = False
# ----------------
if launch_pad_activation:
start_vel = blast_off_start_traj(100, 50)
else:
start_vel = blast_off_start_traj(0, 0)
# A function that plots out the flight trajectory based on the start velocity given.
plot_trajectory(start_vel)
plt.show()