-
Notifications
You must be signed in to change notification settings - Fork 0
/
rocket_flight.py
executable file
·149 lines (123 loc) · 4.58 KB
/
rocket_flight.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#!/usr/bin/env python2
"""
Author: Jonah Miller (jonah.maxwell.miller@gmail.com)
Time-stamp: <2014-09-16 00:28:56 (jonah)>
This is my solution to the rocket flight problem.
Strategy:
We evolve a system of three variables, h, v, and mp, where:
h = the altitude of the rocket
v = the vertical velocity of the rocket
mp = the mass of the propellant remaining on the rocket.
"""
# Imports
# ------------------------------------------------------------------
import numpy
import matplotlib.pyplot as plt
# ------------------------------------------------------------------
# Set the matplotlib rcparams
# -----------------------------------------------------------------
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
# -----------------------------------------------------------------
# Global constants
# ------------------------------------------------------------------
ms = 50 # kg -- Weight of the rocket shell
g = 9.81 # m/s -- Acceleration due to gravity
rho = 1.091 # kg/m^3 -- Average air density
r = 0.5 # m -- Radius of rocket
A = numpy.pi*r**2 # m^2 -- Cross-sectional area of rocket
ve = 325 # m/s -- Exhaust speed
CD = 0.15 # -- Drag coefficient
mp0 = 100 # kg -- Initial weight of the rocket propellant
dt = 0.1 # s -- time step
# Burn rate of the propellant. Not really a constant, but
# conceptually the same.
burn_rate = lambda time: 20 if time < 5 else 0 # kg/s
# Initial data we choose. Assuming initial height and velocity are zero
# but this isn't really necessary.
h0 = 0.
v0 = 0.
# ------------------------------------------------------------------
# Functions
# ------------------------------------------------------------------
# Our right-hand-side function such that du/dt = f(u) for some vector
# u. Our vector u is (h,v,mp).
def rhs(u,t):
"""
Returns the right-hand-side of the differential equation
du/dt = f(u,t)
for the rocket equations of motion
"""
h = u[0] # This is more verboes, but I think it's clearer.
v = u[1]
mp = u[2]
hprime = v # dh/dt
vprime = (-(ms+mp)*g + burn_rate(t)*ve - 0.5*rho*v*abs(v)*A*CD)/(ms+mp)
mpprime = -burn_rate(t)
return numpy.array([hprime,vprime,mpprime])
# We just steal the Euler Step method from the phugoid notebook.
# Only now it depends on time too!
def euler_step(u, t, f, dt):
"""
Returns the solution at the next time-step using Euler's
method.
Parameters
----------
u : array of float
solution at the previous time-step.
t : float
the time at the previous time-step.
f : function
function to compute the right hand-side of the system
of equation.
dt : float
time-increment.
Returns
-------
u_n_plus_1 : array of float
approximate solution at the next time step.
"""
return u + dt * f(u,t)
# A method to perform an Euler solution.
def solve_rocket_equations():
"""
Solves the rocket equations given an initial height h0
and an initial velocity v0.
Uses a first-order forward Euler method with a time step of dt
and evolves from time t = 0 until the rocket crashes
Returns a pair of arrays, the times t and the solutions u(t)
"""
# initial data
u0 = numpy.array([h0,v0,mp0])
t = [0.0] # The array of times
u = [u0] # The array of evolved solutions
# Evolve!
while u[-1][0] >= 0.0:
t.append(t[-1]+dt)
u.append(euler_step(u[-1],t[-1],rhs,dt))
return numpy.array(t),numpy.array(u)
def main():
t,u=solve_rocket_equations()
velocities = u[...,1]
altitudes = u[...,0]
max_velocity = numpy.max(velocities)
mv_index = list(velocities).index(max_velocity)
mv_time = t[mv_index]
mv_h = u[mv_index,0]
max_altitude = numpy.max(altitudes)
ma_index = list(altitudes).index(max_altitude)
ma_time = t[ma_index]
final_time = t[-1]
final_velocity = u[-1,1]
print "We have {} time steps.".format(len(t))
print "At t = 3.2s, the mass in the rocket is: {}".format(u[32,2])
print "The maximum speed of the rocket is: {}".format(max_velocity)
print "This occurs at time t = {} s".format(mv_time)
print "The altitude at this time is h = {} m".format(mv_h)
print "The rocket's maximum altitude is h = {} m".format(max_altitude)
print "This occurs at time t = {} s.".format(ma_time)
print "The rocket hits the ground at t = {} s.".format(final_time)
print "The velocity when it hits the ground is v = {} m/s.".format(final_velocity)
if __name__ == "__main__":
main()