For simplicity: ID
x = v
v = F(x)/m
Errors on Different Time Scales
LTE, GTE
Long-term
Kinetic energy: 1/2ms ||Vs||^2
Potential enery: -G mems/||Xs||
def total_enery(): num_steps = 20000 x = numpy.zeros([num_steps + 1, 2]) v = numpy.zeros([num_steps + 1, 2]) enery = numpy.zeros(num_steps + 1) x[0, 0] = 15e6 x[0, 1] = 1e6 v[0, 0] = 2e3 v[0, 1] = 4e3 for step in range(num_steps): x[step + 1] = x[step] + h * v[step] v[step + 1] = v[step] + h * acceleration(x[step]) return x, energy
Circular Orbit
Orbital period:T
Vs
speed = ||Vs|| = 2πr/T
Acceleration = ||as|| = (2π)^2 r/t^2
Geostationary Orbit
T = 24h
Acceleration => (2π)^2 r/T^2 = 3√GmeT^2/(2π)^2 = 42,000 km
GTE
Error/h1 = E/(h2-h1) = Error = E h1/(h2-h1)
“Order” of a Solves
GTE = C*h order
x = sin(t)
Ch^2new = tolerance
Differential equation
xs = G * mE * -xs / |xs| + G * mm * Xm-xs/ |xm-xs|^3
non-linear
x(t)=sin(t)
x(t)=sin(x(t))
x(t)=cos(y(t))
t=(y(t))^2
differential equation that govern the motion
simple solution method
import math
from xxx import *
moon_distance = 384e6
def orbit():
num_steps = 50
x = numpy.zeros([num_steps + 1, 2])
for i in range(num_steps + 1):
angle = 2.*path.pi * i / num_steps
x[i, 0] = moon_distance * math.cos(angle)
x[i, 1] = moon_distance * math.sin(angle)
return x
x = orbit()
@show_plot
def plot_me():
matplotlib.pyplot.axis('equal')
matplotlib.pyplot.plot(x[:, 0], x[:, 1])
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Longitudinal position in m')
axes.set_ylabel('Lateral position in m')
plot_me()
import math
from xxx import *
h = 0.1
g = 9.81
acceleration = numpy.array
initial_speed = 20. # m / s
@show_plot
def trajectory():
angles = numpy.linspace(20., 70., 6)
num_steps = 30
x = numpy.zeros([num_steps + 1, 2])
v = numpy.zeros([num_steps + 1, 2])
for angle in angles:
matplotlib.pyplot.plot(x[:, 0], x[:, 1])
matplotlib.pyplot.axis('equal')
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Horizontal position in m')
axes.set_ylabel('Vertical position in m')
return x, v
trajectory()
Vector
vector is arrows
a-> = (4 3)
b-> = (2 1)
c-> = (1 -3)
b + c = sum of vector
||a|| = √4^2 + 3^2 = 5
Newton’s low of gravity
F = G m1*m2 / d^2
Gravity = 6.67*10^-11 Nm^2/kg^2
Mass of the moon
m2 = 1.62 m/s^2
R = 1.737*10^6m
A2 = 4A, A1
xM, dES, dMS
import numpy earth_mass = 5.97e24 moon_mass = 7.35e22 gravitational_constant = 6.67e-11 def acceleration(moon_position, spaceship_position): vector_to_moon = moon_position - spaceship_position vector_to_earth = - spaceship_position nonsense = numpy.linalg.norm(moon_position) * spaceship_position return nonsense
Gravity Acceleration
F = m * a(9.81m/s^2)
x(t) = v(t)
v(t) = a(t) = F/m
small time steps of size h
x(h) = x(0) + hv(0)
v(h) = v(0) + h F/m
from xxxplots import
def forward_euler():
h = 0.1
g = 9.81
num_steps = 50
t = numpy.zeros(num_steps + 1)
x = numpy.zeros(num_steps + 1)
v = numpy.zeros(num_steps + 1)
for step in range(num_steps):
t[step + 1] = h * (step + 1)
x[step + 1] = x[step] + h * v[step]
v[step + 1] = v[step] - h * g
return t, x, v
t, x, v = forward_euler()
@show_plot
def plot_me():
axes_height = matplotlib.pyplot.subplot(211)
matplotlib.pyplot.plot(t, x)
axes_velocity = matplotlib.pyplot.subplot(212)
matplotlib.pyplot.plot(t, v)
axes_height.set_ylabel('Height in m')
axes_velocity.set_ylabel('Velocity in m/s')
axes_velocity.set_xlabel('Time in s')
# Uncomment the line below when running locally.
# matplotlib.pyplot.show()
plot_me()
import numpy import matplotlib.pyplot def forward_euler(): h = 0.1 g = 9.81 num_steps = 50 t = numpy.zeros(num_steps + 1) x = numpy.zeros(num_steps + 1) v = numpy.zeros(num_steps + 1) for step in range(num_steps): t[step + 1] = h * (step + 1) x[step + 1] = x[step] + h * v[step] v[step + 1] = v[step] - h * g return t, x, v t, x, v = froward_euler()
if so
k = 130
m = 7655
if k > 42:
print 'Hi'
elif m > 43:
print 'Hello'
else:
print 'Bye'
======================== RESTART: C:/Python27/test.py ======================== Hi >>>
import numpy n = numpy.zeros(5) n[0] = 20. n[4] = 99. print n print n[-1]
import numpy p = numpy.zeros([2, 3]) p[0, 1] = 7. p[1, 2] = 8. print p p = 1. + 2. * p print p
def three_times(u):
r = 2 * u
return r + u
print three_times(1), three_times(7)
======================== RESTART: C:/Python27/test.py ======================== 3 21 >>>
import math
from xxx import *
def sin_cos():
num_points = 50
x = numpy.zeros(num_points)
sin_x = numpy.zeros(num_points)
cos_x = numpy.zeros(num_points)
for i in range(num_points):
x[i] = 2. * math.pi * i / (num_points - 1.)
sin_x[i] = math.sin(x[i])
cos_x[i] = math.cos(x[i])
return x, sin_x, cos_x
x, sin_x, cos_x = sin_cos()
@show_plot
def plot_me():
matplotlib_pyplot.plot(x, sin_x)
matplotlib_pyplot.plot(x, cos_x)
plot_me()
Python Playground
import math f = math.sqrt(4) g = radius_of_earth * math.sin(math.pi / 180.0 * 23.4) g = radius_of_earth * math.sin(math.pi / 180.0 * 23.4) print f, g
======================== RESTART: C:/Python27/test.py ======================== 2.0 2530229.21123 >>>
for h in range(7):
i = h**2
print i
======================== RESTART: C:/Python27/test.py ======================== 0 1 4 9 16 25 36
k = 42
m = 43
while k < 130 or m == 140: k += 1 m += k print k, m [/python]
======================== RESTART: C:/Python27/test.py ======================== 130 7655 >>>
velocity
Velocity = Derivative = Rate of Change
v(t)=dx(t)/dt = x'(t) = x(t)
x(t)=(t+1)^3
(t+h+1)^3 – (t+1)^3 /h
Newton’s Second Law of motion
Force = mass * αacceleration
V = m/s / s = m / s^2
answer = 42 radius_of_earth = 6.371e6 #test number = 43.0 #13test = 44.4 a = 3 * answer + 4 b = radius_of_earth ** 2 + (a + 3) ** 0.5 answer += 7 c = 5 / 4 d = 5 / 4. e = 1.0 * 5 / 4 print a, b, answer, c, d, e
======================== RESTART: C:/Python27/test.py ======================== 130 4.0589641e+13 49 1 1.25 1.25 >>>
Houston We Have a Problem
Apollo 13
TLI, MCC-2, DPS-1
what is important?
mass of the spacecraft, size of the moon, motion of the moon, size of the earth, motion of the earth, 3D
sine & friends
b/c = sinβ
a/c = cosβ
60°=1rad、0.1rad=6°、sin40°=6.4
sin(40)=0.6427876096865393
Motion of the moon around the earth
t = time
4.10^8 2πt*27days