Authors: Sreenidhi Chalimadugu, Rachel Hwang

Views : 125

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)