import plotly.express as px
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

Mission Milestones

This is the 4th notebook in the series. Earlier parts:

• ./010-video-export.html - for exporting still frames and timecodes for the videos
• ./020-video-analysis.html - for analyzing the video frames to match mission time with key events
• ./030-flight-data.html - for analyzing the flight data

Here we will compare the flight data with the mission timeline

df = pd.read_csv("mission-timeline.csv")
df

	Unnamed: 0	Frame	Time (s)	Description	∆T (s)
0	0	25	18.145000	First frame showing rocket plume	0.000000
1	1	31	18.245000	Rocket has cleared the viewport	0.100000
2	2	91	19.245000	Bottle re-acquired	1.100000
3	3	176	20.663333	Rocket at apogee, starting to descend	2.518333
4	4	208	21.196667	Rocket clearly descending and leaning right	3.051667
5	5	245	21.813333	Rocket starts executing flip to nose down	3.668333
6	6	265	22.146667	Rocket pointing downwards, precessing around z...	4.001667
7	7	388	24.198333	Possible rocket crash into building	6.053333

We also found that the launch moment according to the onboard clock was 442461 ms.

LAUNCH_TIME = 442461

Let's look how the video compares with the data

We can load in the sensor data and plot the key events from the video

baro = pd.read_csv('data/BarometerData-20250405-162645.csv')
baro = baro.iloc[11:]
baro = baro[baro['timestamp'] > 440000]
baro['mission_time'] = (baro.timestamp - LAUNCH_TIME) / 1000
baro = baro.set_index('mission_time')
baro.head()

	time	timestamp	pressure	altitude	average_altitude	temperature
mission_time
-2.428	16:20:54.830000	440033	102139.593750	0.210183	0.046781	28.620001
-2.175	16:20:55.083000	440286	102143.968750	-0.150846	0.055706	28.629999
-1.923	16:20:55.335000	440538	102141.726562	0.034192	0.055847	28.650000
-1.669	16:20:55.589000	440792	102138.843750	0.272032	0.047770	28.629999
-1.419	16:20:55.839000	441042	102140.531250	0.132747	0.044345	28.639999

imu = pd.read_csv('data/AccelData-20250405-162645.csv')
imu = imu.iloc[1000:]
imu = imu[imu['timestamp'] > 440000]
imu['mission_time'] = (imu.timestamp - LAUNCH_TIME) / 1000
imu = imu.set_index('mission_time')
imu['total_accel'] = np.sqrt(imu['accelX']**2 + imu['accelY']**2 + imu['accelZ']**2)
imu.head()

	time	timestamp	accelX	accelY	accelZ	gyroX	gyroY	gyroZ	temperature	total_accel
mission_time
-2.459	16:20:54.799000	440002	-0.081984	-0.155672	-1.013576	-1.40	-2.24	2.24	42.7500	1.028733
-2.451	16:20:54.807000	440010	-0.083448	-0.093208	-0.195688	-1.19	-2.59	1.82	42.7500	0.232261
-2.443	16:20:54.815000	440018	-0.122976	-0.087352	-1.002840	-1.82	-2.73	-8.12	42.6875	1.014121
-2.435	16:20:54.823000	440026	-0.131272	-0.179584	-0.770064	-2.80	-0.91	-4.20	42.6875	0.801549
-2.420	16:20:54.838000	440041	0.058072	-0.064904	-1.726544	0.49	-1.96	12.88	42.6875	1.728739

Let's create a function to annotate the plot with the mission milestones

def annotate_plot(fig, y):

    for _, row in df.iterrows():
        fig.add_vline(
            x=row['∆T (s)'],
            line_dash="dash",
            line_color="red",
        )
        fig.add_annotation(
            x=row['∆T (s)'],
            y=y,
            text=row['Description'],
            showarrow=False,
            textangle=-90,
            font=dict(color="grey", size=10),
            xanchor="left",
            yanchor="top"
        )

Flight profile

Let's first examine how our "video" milestones compare to the real flight data

filtered_data = baro['altitude']

fig = px.scatter(x=filtered_data.index, y=filtered_data.values, 
                 title="Rocket Altitude Measurements",
                 labels={"x": "Time (s)", "y": "Altitude (m)"})

fig.update_traces(mode='lines+markers')

fig.update_layout(height=1000)
fig.update_xaxes(range=[-0.2, 7])
fig.update_yaxes(range=[-1, 100])

annotate_plot(fig, 90)

fig.show()

Ok - so we weren't particularly good at estimating the apogee, but we were close. It also from this data shows that the rocket didn't crash into the building, and we have data all the way to the ground, although there is possibly a slightly slowdown in descent for the last few measurements which COULD indicate we hit a tree or something, but not clear from these data.

The first "in flight" datapoint is estimating higher altitude than what seems to fit the curve. Apart from being under quite a lot of stress at the time, the bottle has just "deflated" which could cause a bit of underpressure in the "payload fairing" which sits on top of the bottle. It is reasonable to assume that this point is unreliable and that the actual altitude is lower than indicated.

Accelerometer Data

This should be the clearest dataset for the "BIG" mission milestones: launch and crash. Let's see how it lines up.

filtered_data = imu['total_accel']

fig = px.scatter(x=filtered_data.index, y=filtered_data.values, 
                 title="Rocket Acceleration",
                 labels={"x": "Time (s)", "y": "Acceleration (G)"})

fig.update_traces(mode='lines')

fig.update_layout(height=1000)
fig.update_xaxes(range=[-0.2, 7])
fig.update_yaxes(range=[-1, 25])

annotate_plot(fig, 24)

fig.show()

Prior to launch we see that accelleration is about 1G, which is what we expect.

After launch, there is a massive spike for about 0.1 second, which presumably corresponds to the water ballast being ejected.

After this, the rocket has a falling thrust from about 1G to 0G is the air in the bottle equalizes with the outside.

As we observe the 'tumble' of the rocket as it transitions to nose-down we see a slightly bump in the acceleration, which must correspond to increased drag during the tumble.

Finally, at about 6.5s of flight, we see the rocket hit the ground, possibly bounce, and then perhaps roll.

We suffer a loss of signal at 6.95s caused by a reset of the flight computer.

So far we have looked at total acceleration, but we can also look at the individual components of acceleration.

filtered_data = imu[["accelX", "accelY", "accelZ"]].copy()
filtered_data["Time"] = filtered_data.index  # Turn index into a column for plotly

# Melt the DataFrame to long format
long_data = filtered_data.melt(id_vars="Time", var_name="Axis", value_name="Acceleration")

fig = px.line(
    long_data,
    x="Time",
    y="Acceleration",
    color="Axis",
    title="Rocket Acceleration - Each Axis",
    labels={"Time": "Time (s)", "Acceleration": "Acceleration (G)"}
)

fig.update_traces(mode='lines')

fig.update_layout(height=1000)
fig.update_xaxes(range=[-0.2, 7])
fig.update_yaxes(range=[-16.5, 16.5])

annotate_plot(fig, 15)

fig.show()

A few key observations here:

1. Acceleration due to gravity is about 1G as seen prior to launch - so this means our instruments are properly calibrated.
2. During launch, the acceleration in the Z axis (along the rocket's axis) is clipped at -16G which is the limit of the IMU. Launch acceleration has clearly exceeded this limit.
3. Post launch we see positive z-axis acelleration which is most likely due to aerodynamic drag on the rocket. This matches the observations in that the drag is reducing as the rocket slows down, then raises again as the rocket speeds up. We see a bump as the rocket flips, and the drag is always towards the rear of the rocket (+z axis) corresponding to a retarding force.

filtered_data = imu[["accelX", "accelY", "accelZ"]].copy()
filtered_data["Time"] = filtered_data.index  # Turn index into a column for plotly

# Melt the DataFrame to long format
long_data = filtered_data.melt(id_vars="Time", var_name="Axis", value_name="Acceleration")

fig = px.line(
    long_data,
    x="Time",
    y="Acceleration",
    color="Axis",
    title="Rocket Acceleration - Each Axis",
    labels={"Time": "Time (s)", "Acceleration": "Acceleration (G)"}
)

fig.update_traces(mode='lines')

fig.update_layout(height=1000)
fig.update_xaxes(range=[-0.2, 7])
fig.update_yaxes(range=[-0.5, 1.5])

annotate_plot(fig, 1.4)

fig.show()

So, what is immediately obvious is that the initial retarding acceleration is higher than during the descent phase. This indicates that the launch velocity was higher than achieved during the descent phase.

Unexpectedly the rocket is also seeing some backward acceleration close to apogee. It could be that the rocket has started falling "ass first" at this stage, before flipping, but it's also worth noting that the magnitude of these forces are similar to what is seen along the x and y axes caused by slight tumbling of the rocket.

It could indicate that the rocket reaches apogee earlier than we estimated (an in fact more consistent with the video data than the barometric data), and one could postulate that there is possibly a slight delay in the barometric data due to over/underpressure within the payload fairing. This hypothesis is also supported by the measurements of the barometric altitude contiuing to fall even after we know the rocket has impacted the ground.

MECO to apogee

Eyeballing the plot above, it seems that the curve is not linear, but rather follows an exponential decay. This is typical for objects subjected to drag forces at high velocities, or when the air is turbulent.

• Shape and General Trend: The data show a clear initial peak followed by a rapid decline that gradually tapers off. This shape is typical for objects subjected to drag forces, which are usually proportional to the square of the velocity (quadratic drag).

• Suggested Model: Based on visual inspection, the data likely follows a quadratic drag model. Air resistance (drag) at higher velocities is commonly modeled as:

$$ F_{\text{drag}} = \frac{1}{2}\rho C_d A v^2 $$

where:

• $ \rho $ is air density,

• $ C_d $ is the drag coefficient,

• $ A $ is the cross-sectional area,

• $ v $ is the velocity.

• Curve Fit Expectation: The data should be fit well by a model resembling:

$$ y(t) = a e^{-b t} + c $$

or, more directly tied to velocity-based drag force:

$$ F_{\text{drag}}(t) \propto v(t)^2 $$

from scipy.optimize import curve_fit

filtered_data = imu[["accelZ"]].copy()

# All data in launch phase and a bit beyond
filtered_data = filtered_data[(filtered_data.index > 0) & (filtered_data.index < 3)] 

# Only forces pushing "down" 
filtered_data = filtered_data[filtered_data["accelZ"] >= 0] 

# We adjust time to be zero at the time of MECO when the rocket starts it's ballistic phase
t_start = min(filtered_data.index)
filtered_data.index = filtered_data.index - t_start


# Exponential decay model function
def model_func(t, a, b, c):
    return a * np.exp(-b * t) + c

# Extracting data for fitting
time = filtered_data.index
force = filtered_data['accelZ'].values

# Initial guess for the parameters (a, b, c)
initial_guess = [force.max(), 1.0, force.min()]

# Curve fitting
params, covariance = curve_fit(model_func, time, force, p0=initial_guess)

# Fitted parameters
a_fit, b_fit, c_fit = params
print(f"Fitted Parameters:\na = {a_fit}\nb = {b_fit}\nc = {c_fit}")
print(f"Max time: {max(time):.3f}s")

# Plotting the original data and fitted curve
plt.figure(figsize=(10, 6))
plt.scatter(time, force, label='Measured Data', s=10, color='blue')
plt.plot(time, model_func(time, *params), label='Fitted Curve', color='red', linewidth=2)
plt.xlabel('Time')
plt.ylabel('Acceleration (G)')
plt.title('Curve Fit of Acceleration (Retardation) vs. Time')
plt.legend()
plt.grid()

Fitted Parameters:
a = 1.317398309564944
b = 1.0491244887060678
c = -0.10783706262291398
Max time: 2.376s

The total forces acting on the rocket are now given by the curve above plus gravity. At apogee the rocket's velocity has fallen to zero.

$$ a_{\text{drag}}(t_{\text{stop}}) = 0 \Rightarrow a e^{-b t_{\text{stop}}} + c = 0 $$

We can calculat the time of v=0:

$$ t_{\text{stop}} = -\frac{1}{b} \ln\left(-\frac{c}{a}\right) $$

The drag acceleration causes deceleration, so:

$$ \frac{dv}{dt} = -g - a_{\text{drag}}(t) = -g - a e^{-bt} - c $$

Then:

$$ v_0 = \int_0^{t_{\text{stop}}} (g + a e^{-bt} + c) \, dt $$

In units of G, ( g = 1 ), so:

$$ v_0 = \int_0^{t_{\text{stop}}} (1 + a e^{-bt} + c) \, dt = (1 + c)t_{\text{stop}} + \frac{a}{b}(1 - e^{-b t_{\text{stop}}}) $$

# Fitted parameters (from your curve fit)
a = 1.317398309564944
b = 1.0491244887060678
c = -0.10783706262291398

# Step 1: Find time when a_drag(t) = 1 G
t_stop = -np.log(-c / a) / b

# Step 2: Compute v0 in G·s
v0_g = (1 + c) * t_stop + (a / b) * (1 - np.exp(-b * t_stop))

# Convert to m/s (optional)
v0_ms = v0_g * 9.81
v0_kmh = v0_ms * 3.6

# Step 3: Compute altitude h(t_stop)
h_stop = ((1 + c) / 2) * t_stop**2 + (a / b) * t_stop + (a / b**2) * (np.exp(-b * t_stop) - 1)

# Convert to meters (from G * s²)
h_stop_m = h_stop * 9.81  # because 1 G = 9.81 m/s²

print(f"Time when velocity = 0: {t_stop:.4f} s (or {t_stop+t_start:.4f} in mission time)")
print(f"Initial velocity: {v0_g:.4f} G·s or {v0_ms:.2f} m/s or {v0_kmh:.2f} km/h")
print(f"Altitude at v = 0: {h_stop:.4f} G·s² or {h_stop_m:.2f} m")

Time when velocity = 0: 2.3856 s (or 2.4996 in mission time)
Initial velocity: 3.2813 G·s or 32.19 m/s or 115.88 km/h
Altitude at v = 0: 4.4354 G·s² or 43.51 m

From analysis of video data in ./060-launch-analysis.html we estimated launch velocity to be 35.56 m/s 🎉

Apogee to "landing"

Now, let's see what happens after apogee. It seems the rocket initially starts falling backwards, but then flips over and starts to fall nose first.

Because the rocket flips during this phase, it makes more sense to only consider the magnitude of the acceleration.

from scipy.optimize import curve_fit

filtered_data = imu[["accelZ"]].copy()

# All data in launch phase and a bit beyond
filtered_data = filtered_data[(filtered_data.index > 2.4996) & (filtered_data.index < 6.515)] 

filtered_data['accelZ'] = filtered_data['accelZ'].abs()
# filtered_data = filtered_data[(filtered_data.index < 3.4) | (filtered_data.index > 4.5)]
filtered_data.plot()

# display(filtered_data.index)

<Axes: xlabel='mission_time'>

Let's snip out the bits around the flip - let's say from 3.5 to 4.5 seconds.

# Delete rows where index is between 3.3 and 4.6
filtered_data = filtered_data[(filtered_data.index < 3.3) | (filtered_data.index > 4.6)]
filtered_data.plot(style='.')

<Axes: xlabel='mission_time'>

t_start = min(filtered_data.index)
filtered_data.index = filtered_data.index - t_start

# Extracting data for fitting
time = filtered_data.index
force = filtered_data['accelZ'].values

def model_func(t, a, b):
    return a * (np.exp(b * t) - 1)

# Curve fitting with the adjusted model
initial_guess = [force.max(), 1.0]

params, covariance = curve_fit(model_func, time, force, p0=initial_guess)

# # Fitted parameters
a_fit, b_fit = params
print(f"Fitted Parameters:\na = {a_fit}\nb = {b_fit}")
print(f"Max time: {max(time):.3f}s")

# # Plotting the original data and fitted curve
plt.figure(figsize=(10, 6))
plt.scatter(time, force, label='Measured Data', s=10, color='blue')
plt.plot(time, model_func(time, *params), label='Fitted Curve', color='red', linewidth=2)
plt.xlabel('Time')
plt.ylabel('Acceleration (G)')
plt.title('Curve Fit of Acceleration vs. Time')
plt.legend()
plt.grid()

Fitted Parameters:
a = 0.1312788063747338
b = 0.4175709077175023
Max time: 4.000s

import numpy as np
from scipy.integrate import cumulative_trapezoid
import matplotlib.pyplot as plt

# Constants
g = 9.81  # gravitational acceleration (m/s^2)

# Define time range
t = np.linspace(0, max(time), 1000)

# Acceleration due to air resistance
def air_resistance_accel(t):
    return a_fit * (np.exp(b_fit * t) - 1)

# Total acceleration
def total_accel(t):
    return g - air_resistance_accel(t)

accel_total_values = total_accel(t)

# Velocity and displacement
velocity = cumulative_trapezoid(accel_total_values, t, initial=0)
displacement = cumulative_trapezoid(velocity, t, initial=0)

# Plotting
fig, axes = plt.subplots(3, 1, figsize=(10, 15))

# Acceleration plot
axes[0].plot(t, accel_total_values, label='Total Acceleration')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Acceleration (m/s²)')
axes[0].set_title('Total Acceleration vs Time')
axes[0].grid(True)
axes[0].legend()

# Velocity plot
axes[1].plot(t, velocity, color='orange', label='Velocity')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Velocity (m/s)')
axes[1].set_title('Velocity vs Time')
axes[1].grid(True)
axes[1].legend()

# Displacement plot
axes[2].plot(t, displacement, color='green', label='Displacement')
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Displacement (m)')
axes[2].set_title('Displacement vs Time')
axes[2].grid(True)
axes[2].legend()

plt.tight_layout()
plt.show()

Gyroscope Data

Let's look at the gyroscope data to see if we can get any more information about the tumble.

filtered_data = imu[["gyroX", "gyroY", "gyroZ"]].copy()
filtered_data["Time"] = filtered_data.index  # Turn index into a column for plotly

# Melt the DataFrame to long format
long_data = filtered_data.melt(id_vars="Time", var_name="Axis", value_name="Gyro")

fig = px.line(
    long_data,
    x="Time",
    y="Gyro",
    color="Axis",
    title="Rocket Gyroscope Data",
    labels={"Time": "Time (s)", "Gyro": "Gyroscope Value (deg/s)"}
)
fig.update_layout(height=1000)
fig.update_xaxes(range=[-0.2, 7])
fig.update_yaxes(range=[-800, 800])

annotate_plot(fig, 800)

fig.show()

Interestingly during the thrust phase, the rocket is spinning slightly in the positive direction, but as soon as the main thrust drops off, the fins start to induce a strong rotation which also helps dampen any pitch or yaw oscillations. As the rocket slows down towards apogee, the rotation slows down, but it maintains rotation through the slowest part of the flight.

Att the tumble point, we see both the pitch and yaw rotation increase, which is what we expect, and they settle into a slightly oscillation which is also dampened as the velocity of the rocket increases during the descent phase and the fins again increase roll rotation.

Finally, as the rocket hits the ground, we see a large, chaotic spike in the rotation, which is likely due to the impact.