Navigation
Lab 10
Information about the lab goes here.
from __future__ import division
from visual import *
"""
Filename: YourName_Lab10_PHYS495_1.py
Author: Your Name
Date: 04/18/2007
Description: A lunar mission simulation.
We use SI units.
Method: Velocity-Verlet.
"""
#---------------------------------------------------------
# Useful quantities
#---------------------------------------------------------
Earth_radius = 6378.1e3 # meters
Earth_mass = 5.9736e24 # kilograms
Moon_radius = 1.7371e6 # meters
Moon_mass = 7.349e22 # kilograms
Moon_orbital_radius = 384.4e6 # meters
Moon_orbital_velocity = 2*pi*Moon_orbital_radius / (27.3 * 86400) # m/s
#*********** PARAMETERS TO EDIT *********************
#
# Change the following parameters to get your
# rocket to land on the Moon.
#
# Moon_angle
#
# Thrust_angle (change the number added to Latitude)
#
# Angle_rotation_rate (change the last number--in the
# denominator)
#
#****************************************************
Moon_angle = 0.0 # degrees; Moon's initial position
Latitude = 0 #28.5 # degrees--- DON'T CHANGE THIS
Thrust_angle = float(Latitude + 8) # degrees
Angle_rotation_rate = (90 - Thrust_angle) / 90 # degrees / second
G = 6.673e-11 # SI units
t = 0
t_max = 7 * 86400 # seconds
dt = 1.0 # seconds
dt_2 = dt / 2 # half step
##t_max = 2 * dt
pi = acos(-1.0)
deg_2_rad = pi / 180
#---------------------------------------------------------
# Scene attributes
#---------------------------------------------------------
w = 800
h = 600
scene.title = "VPython: MOON MISSION 1 (Lab 10)"
scene.x = 0
scene.y = 0
scene.width = w
scene.height = h
scene.range = (62 * Earth_radius, 62 * Earth_radius, 62 * Earth_radius)
scene.autoscale = 0 ##0 means autoscaling is OFF
scene.userzoom = 1 ##0 means user cannot zoom
scene.userspin = 1 ##0 means user cannot spin
scene.lights = [vector(0,0,1)]
scene.ambient = 0.5
def check_for_pause(scene): #checks for pause request
if scene.mouse.clicked:
scene.mouse.getclick()
pause(scene)
def pause(scene):
while 1:
if scene.mouse.clicked:
scene.mouse.getclick()
break
data_scene = display()
data_scene.title = "Flight Data (Lab 10)"
data_scene.width = w
data_scene.height = h
data_scene.x = scene.x + 805
data_scene.y = scene.y
scene.select()
#---------------------------------------------------------
# Create the objects for the animation
#---------------------------------------------------------
earth = sphere()
earth.mass = float(Earth_mass)
earth.radius = float(Earth_radius)
earth.pos = vector(float(0), float(0), float(0))
earth.vel = vector(float(0), float(0), float(0))
earth.color = color.green
equator = ring()
equator.pos = vector(float(0), float(0), float(0))
equator.axis = vector(0,0,1)
equator.radius = 1.00 * earth.radius
equator.thickness = 0.01 * earth.radius
moon = sphere()
moon.mass = float(Moon_mass)
moon.radius = float(Moon_radius)
moon.pos = vector(Moon_orbital_radius * cos(Moon_angle * deg_2_rad),
Moon_orbital_radius * sin(Moon_angle * deg_2_rad),
0.0)
moon.vel = (Moon_orbital_velocity *
vector(-1*sin(Moon_angle * deg_2_rad),
cos(Moon_angle * deg_2_rad),
0.0))
moon.force = (-1 * G * (moon.pos - earth.pos) * moon.mass * earth.mass /
mag(moon.pos - earth.pos)**3)
moon.acc = moon.force / moon.mass
moon.color = color.white
moon.orbit = curve(color = color.green)
rocket = sphere()
rocket.orbit = curve(color = color.red)
rocket.color = color.blue
rocket.radius = 0.05 * earth.radius
rocket.mass = 2.659e6 #kilograms
rocket.original_mass = 2.659e6 #kilograms
rocket.stage_1_mass = 2.08e6 # kilograms
rocket.stage_2_mass = 4.28e5 # kilograms
rocket.stage_3_mass = 1.05e5 # kilograms
rocket.stage_1_dmdt = 13000 # 12988.505747 # kg/s
rocket.stage_2_dmdt = 1070 # 1231.16883117 # kg/s
rocket.stage_3_dmdt = 210 # 245.901639344 # kg/s
rocket.stage_1_v = 2615.38461538 # m/s
rocket.stage_2_v = 4672.8971926 # 4065 # m/s
rocket.stage_3_v = 4761.9047619 # 4065 # m/s
rocket.stage_1_thrust = rocket.stage_1_dmdt * rocket.stage_1_v *\
vector(float(cos(Thrust_angle * deg_2_rad)),
float(sin(Thrust_angle * deg_2_rad)),
float(0.0)) # newtons
rocket.stage_2_thrust = vector(0.0, 0.0, 0.0) # newtons
rocket.stage_3_thrust = vector(0.0, 0.0, 0.0) # newtons
rocket.pos = vector(float(earth.radius * cos(Latitude * deg_2_rad)),
float(0.0),
float(earth.radius * sin(Latitude * deg_2_rad)))
rocket.vel = vector(float(0.0),
float(2*pi*Earth_radius*cos(Latitude)/86400),
float(0.0))
rocket.weight = (-1 * G * (rocket.pos - earth.pos) * rocket.mass * earth.mass
/mag(rocket.pos - earth.pos)**3
+ -1 * G * (rocket.pos - moon.pos) * rocket.mass * moon.mass
/mag(rocket.pos - moon.pos)**3)
rocket.acc = (rocket.weight + rocket.stage_1_thrust) / rocket.mass
#---------------------------------------------------------
# Create a label for data
#---------------------------------------------------------
orbital_data = label()
orbital_data.display = data_scene
orbital_data.pos = vector(0, 0, 0)
orbital_data.xoffset = 0
orbital_data.yoffset = -110
orbital_data.line = 0
orbital_data.opacity = 0
orbital_data.box = 0
orbital_data.text = ("H (km) = %3.6f \n"
"To go (km) = %3.6f \n"
"v = %3.6f \n"
"a = %3.6f \n"
"Thrust angle = %3.3f\n"
"stage 1 mass = %3.6e\n"
"stage 2 mass = %3.6e\n"
"stage 3 mass = %3.6e\n"
"total mass = %3.6e\n"
"total mass (per of original) = %3.3f\n"
"stage 1 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"stage 2 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"stage 3 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"weight = %3.6e\n"
"t = %3.3f \n"
"t (hr) = %3.3f\n"
%(((mag(rocket.pos) - earth.radius)/1000,
mag(rocket.pos - moon.pos)/1000,
mag(rocket.vel),
mag(rocket.acc),
Thrust_angle,
rocket.stage_1_mass,
rocket.stage_2_mass,
rocket.stage_3_mass,
rocket.mass,
100 * rocket.mass / rocket.original_mass,
rocket.stage_1_thrust.x,
rocket.stage_1_thrust.y,
rocket.stage_1_thrust.z,
rocket.stage_2_thrust.x,
rocket.stage_2_thrust.y,
rocket.stage_2_thrust.z,
rocket.stage_3_thrust.x,
rocket.stage_3_thrust.y,
rocket.stage_3_thrust.z,
mag(rocket.weight),
t,
t/3600)))
#---------------------------------------------------------
# The animation
#---------------------------------------------------------
"""
The animation loop begins with a pause statement.
This gives the user time to find the animation window
on the Windows desktop before the animation begins.
The user needs to click the animation scene to start
the animation.
"""
pause(scene)
while(abs(t_max - t) > dt_2):
## rate(300) # increase number to increase speed of the animation
check_for_pause(scene)
t += dt # increment the elapsed time
## if t > 14000: # Recenter the scene to follow the rocket
## scene.center = rocket.pos
if mag(rocket.pos - moon.pos) < moon.radius:
t_max = t
moon.vel += moon.acc * dt_2
moon.pos += moon.vel * dt
moon_force = ( -1 * G * (moon.pos - earth.pos) * moon.mass * earth.mass /
mag(moon.pos - earth.pos)**3)
moon.orbit.append(pos = moon.pos)
moon.acc = moon_force / moon.mass
moon.vel += moon.acc * dt_2
rocket.vel += rocket.acc * dt_2
rocket.pos += rocket.vel * dt
rocket.orbit.append(pos = rocket.pos)
rocket.weight = (-1 * G * (rocket.pos - earth.pos) * rocket.mass * earth.mass
/mag(rocket.pos - earth.pos)**3
+ -1 * G * (rocket.pos - moon.pos) * rocket.mass * moon.mass
/mag(rocket.pos - moon.pos)**3)
if t <= 160 : # Stage 1 (160 sec)
rocket.stage_1_mass -= rocket.stage_1_dmdt * dt
rocket.mass -= rocket.stage_1_dmdt * dt
rocket.stage_1_thrust = rocket.stage_1_dmdt * rocket.stage_1_v *\
vector(float(cos(Thrust_angle * deg_2_rad)),
float(sin(Thrust_angle * deg_2_rad)),
float(0.0)) # newtons
rocket.acc = (rocket.stage_1_thrust + rocket.weight) / rocket.mass
elif 160 < t <= 560 : # Stage 2 (400 sec)
if Thrust_angle <= 90.0:
Thrust_angle += Angle_rotation_rate
rocket.stage_1_thrust = vector(0.0, 0.0, 0.0)
rocket.stage_2_mass -= rocket.stage_2_dmdt * dt
rocket.mass -= rocket.stage_2_dmdt * dt
rocket.stage_2_thrust = rocket.stage_2_dmdt * rocket.stage_2_v *\
vector(float(cos(Thrust_angle * deg_2_rad)),
float(sin(Thrust_angle * deg_2_rad)),
float(0.0)) # newtons
rocket.acc = (rocket.stage_2_thrust + rocket.weight) / rocket.mass
elif 560 < t <= 700 : # Stage 3, first burn (140 sec)
rocket.stage_2_thrust = vector(0.0, 0.0, 0.0)
rocket.stage_3_mass -= rocket.stage_3_dmdt * dt
rocket.mass -= rocket.stage_3_dmdt * dt
rocket.stage_3_thrust = rocket.stage_3_dmdt * rocket.stage_3_v *\
(vector(rocket.vel)/mag(rocket.vel)) # newtons
rocket.acc = (rocket.stage_3_thrust + rocket.weight) / rocket.mass
elif 700 < t <= 13500 : # Coast period
rocket.stage_3_thrust = vector(0.0, 0.0, 0.0)
rocket.acc = (rocket.stage_3_thrust + rocket.weight) / rocket.mass
elif 13500 < t <= 13860: #13860 : # Second and final third stage burn (360 sec)
rocket.stage_3_mass -= rocket.stage_3_dmdt * dt
rocket.mass -= rocket.stage_3_dmdt * dt
rocket.stage_3_thrust = rocket.stage_3_dmdt * rocket.stage_3_v *\
(vector(rocket.vel)/mag(rocket.vel)) # newtons
rocket.acc = (rocket.stage_3_thrust + rocket.weight) / rocket.mass
else:
rocket.stage_3_thrust = vector(0.0, 0.0, 0.0)
rocket.acc = rocket.weight / rocket.mass
rocket.vel += rocket.acc * dt_2
orbital_data.text = ("H (km) = %3.6f \n"
"To go (km) = %3.6f\n"
"v = %3.6f \n"
"a = %3.6f \n"
"Thrust angle = %3.3f\n"
"stage 1 mass = %3.6e\n"
"stage 2 mass = %3.6e\n"
"stage 3 mass = %3.6e\n"
"total mass = %3.6e\n"
"total mass (per of original) = %3.3f\n"
"stage 1 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"stage 2 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"stage 3 thrust x: %3.6e y: %3.6e z: %3.6e\n"
"weight = %3.6e\n"
"t = %3.3f \n"
"t (hr) = %3.3f\n"
%(((mag(rocket.pos) - earth.radius)/1000,
mag(rocket.pos - moon.pos)/1000,
mag(rocket.vel),
mag(rocket.acc),
Thrust_angle,
rocket.stage_1_mass,
rocket.stage_2_mass,
rocket.stage_3_mass,
rocket.mass,
100 * rocket.mass / rocket.original_mass,
rocket.stage_1_thrust.x,
rocket.stage_1_thrust.y,
rocket.stage_1_thrust.z,
rocket.stage_2_thrust.x,
rocket.stage_2_thrust.y,
rocket.stage_2_thrust.z,
rocket.stage_3_thrust.x,
rocket.stage_3_thrust.y,
rocket.stage_3_thrust.z,
mag(rocket.weight),
t,
t/3600)))