SharedOrbital Mechanics Mini Project.ipynbOpen in CoCalc
Download, Raw, Embed
Authors: Sreenidhi Chalimadugu, Rachel Hwang
Views : 24
In [62]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *
import math
In [63]:
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
kilogram
In [2]:
params = Params(G_const = 6.67408e-11 * m**3/(kg*s**2), m_e = 5.972e24 * kg, m_s = 1.989e30 * kg) # system
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-2-fa5e6484cfa8> in <module> ----> 1 params = Params(G_const = 6.67408e-11 * m**3/(kg*s**2), m_e = 5.972e24 * kg, m_s = 1.989e30 * kg) # system NameError: name 'Params' is not defined
In [65]:
def make_system(params):
    """Makes a System object for the given conditions.
    
    params: Params object
    
    returns: System object
    """
    G_const, m_e, m_s = params.G_const, params.m_e, params.m_s
    
    r_init = 149.6e9 * m
    v_init = 0 *m/s
    init = State(r=r_init, v=v_init)
    t_end = 3600 * 24 * 100 * s
    dt = t_end / 100
    
    return System(params, G_const=G_const, init=init, t_end=t_end, dt=dt)
In [66]:
system = make_system(params)
values
G_const 6.67408e-11 meter ** 3 / kilogram / second ** 2
m_e 5.972e+24 kilogram
m_s 1.989e+30 kilogram
init r 149600000000.0 meter v 0.0 meter / s...
t_end 8640000 second
dt 86400.0 second
In [67]:
def slope_func(state, t, system):
    """Compute derivatives of the state.
    
    state: position, velocity
    t: time
    system: System object
    
    returns: derivatives of y and v
    """
    r, v = state
    G_const, m_s = system.G_const, system.m_s
    
    accel = (G_const * m_s) / r**2
    
    drdt = v
    dvdt = - G_const * m_s / r**2
    
    return drdt, dvdt
In [68]:
slope_func(system.init, 0, system)
(0.0 <Unit('meter / second')>, -0.00593147909577054 <Unit('meter / second ** 2')>)
In [69]:
def event_func(state, t, system):
    r, v = state
    return r
In [70]:
results, details = run_ode_solver(system, slope_func, events=event_func)
details
values
success True
message A termination event occurred.
In [71]:
results.head()
r v
0.0 149600000000.0 meter 0.0 meter / second
86400.0 149577860872.90463 meter -512.4797938745746 meter / second
172800.0 149511436937.47714 meter -1025.263098836451 meter / second
259200.0 149400688842.6613 meter -1538.6543818132357 meter / second
345600.0 149245550867.2312 meter -2052.960096377744 meter / second
In [72]:
results.tail()
r v
5.270400e+06 36542038776.296616 meter -74090.6292983068 meter / second
5.356800e+06 29769553538.081444 meter -84496.11754051299 meter / second
5.443200e+06 21910002441.51194 meter -101602.72652695786 meter / second
5.529600e+06 12099386348.186493 meter -140936.47367582278 meter / second
5.596778e+06 0.0 meter -578088.5928502546 meter / second
In [75]:
def plot_position(results):
    indexes = results.index / (60 * 60 * 24)
    radii = results.r / 1e9
    plot(indexes, radii)
    decorate(xlabel='Time (Days)',
             ylabel='Position (1 billion meters)')
    
plot_position(results)
In [74]:
get_last_label(results) / (3600 * 24)
64.77752403749415

PART 2

In [113]:
def make_system2(params):
    """Makes a System object for the given conditions.
    
    params: Params object
    
    returns: System object
    """
    G_const, m_e, m_s = params.G_const, params.m_e, params.m_s
    
    r_init = Vector(0, 149.6e9) * m
    v_init = Vector(30000, 0) * m/s
    init = State(r=r_init, v=v_init)
    t_end = 3600 * 24 * 365 * s
    dt = t_end / 100
    
    return System(params, G_const=G_const, init=init, t_end=t_end, dt=dt)
In [114]:
system = make_system2(params)
values
G_const 6.67408e-11 meter ** 3 / kilogram / second ** 2
m_e 5.972e+24 kilogram
m_s 1.989e+30 kilogram
init r [0.0 meter, 149600000000.0 met...
t_end 31536000 second
dt 315360.0 second
In [115]:
def slope_func2(state, t, system):
    """Compute derivatives of the state.
    
    state: position, velocity
    t: time
    system: System object
    
    returns: derivatives of y and v
    """
    r, v = state
    G_const, m_s = system.G_const, system.m_s
        
    drdt = v
    dvdt = G_const * -r.hat() * m_s / r.mag**2
    return drdt, dvdt
In [116]:
slope_func2(system.init, 0, system)
(array([30000., 0.]) <Unit('meter / second')>, array([-0. , -0.00593148]) <Unit('meter / second ** 2')>)
In [117]:
def event_func2(state, t, system):
    r, v = state
    return r.mag
In [118]:
results, details = run_ode_solver(system, slope_func2, events=event_func2)
details
values
success True
message The solver successfully reached the end of the...
In [119]:
def plot_position2(results):
    indexes = results.index / (60 * 60 * 24)
    x = results.r.extract('x') / 1e9
    y = results.r.extract('y') / 1e9
    plot(indexes, x)
    plot(indexes, y)

    decorate(xlabel='Time (Days)',
             ylabel='Position (1 billion meters)')
    
plot_position2(results)
In [121]:
def plot_orbit(results):
    indexes = results.index / (60 * 60 * 24)
    x = results.r.extract('x') / 1e9
    y = results.r.extract('y') / 1e9
    plot(x, y)

    decorate(xlabel='Time (Days)',
             ylabel='Position (1 billion meters)')
    
plot_orbit(results)