AAE 532 - Hohmann Transfer gone bad

hw09.py

# AAE 532 - HW09 - Q02
# Author : Thomas Antony ( tantony (at) purdue (dot) edu )
# Code re-purposed from from Udacity.com CS222 - HW02-Q02
# Original code author : Miriam "Swords" Kalk ?? ( I think ) 
# The simulation starts with the spacecraft in an earth-centered
# circular orbit at 175 km altitude, just after finishing a burn which sends it on a 
# Hohmann trajectory to the moon. 
#
# The purpose of the simulation was to examine what happens to the orbit if a capture
# burn was not performed after reaching the moon
#
# Requires : udacityplots.py from https://gist.github.com/3612552, which in turn requires
#            numpy, matplotlib and pyplot libraries
#

import math
from udacityplots import *

earth_mass = 5.97e24 # kg
earth_radius = 6.378e6 # m (at equator)
gravitational_constant = 6.67e-11 # m3 / kg s2
moon_mass = 7.35e22 # kg
moon_radius = 1.74e6 # m
moon_distance = 400.5e6 # m (actually, not at all a constant)
moon_period = 27.3 * 24.0 * 3600. # s
moon_initial_angle = math.pi / 180.  *  24.72611# angle from horizontal

total_duration = 20. * 24. * 3600. # s
marker_time = 0.5 * 3600. # s
tolerance = 100000. # m



def moon_position(time):
# Task 1: Compute the moon's position (a vector) at time t. Let it start at moon_initial_angle, not on the horizontal axis.   
    ###Your code here.
    angle = moon_initial_angle+time*(2*math.pi)/moon_period;
    return numpy.array([moon_distance*math.cos(angle),moon_distance*math.sin(angle)]);

def acceleration(time, position):
# Task 2: Compute the spacecraft's acceleration due to gravity
    moon_pos = moon_position(time);

    mu_e = gravitational_constant*earth_mass;
    mu_m = gravitational_constant*moon_mass;

    moon_dist_vec = position-moon_pos;
    dist_e = numpy.linalg.norm(position)
    dist_m = numpy.linalg.norm(moon_dist_vec)

    acc = -(mu_e / numpy.linalg.norm(position)**3 * position 
                + mu_m / numpy.linalg.norm(moon_dist_vec)**3 * moon_dist_vec);

    return acc  

# Task 5: (First see the other tasks below.) What is the appropriate boost to apply?
# Try -10 m/s, 0 m/s, 10 m/s, 50 m/s and 100 m/s and leave the correct amount in as you submit the solution.

#axes = matplotlib.pyplot.gca()
#axes.set_xlabel('Longitudinal position in m')
#axes.set_ylabel('Lateral position in m')

def apply_boost():

    # Do not worry about the arrays position_list, velocity_list, and times_list.  
    # They are simply used for plotting and evaluating your code, so none of the 
    # code that you add should involve them.
    
    
    boost = 10. # m/s Change this to the correct value from the list above after everything else is done.
    position_list = [numpy.array([0, -6553.14e+3])] # m
    velocity_list = [numpy.array([10936.750, 0])] # m / s
    times_list = [0]
    position = position_list[0]
    velocity = velocity_list[0]
    current_time = 0.
    h = 0.1 # s, set as initial step size right now but will store current step size
    h_new = h # s, will store the adaptive step size of the next step
    mcc2_burn_done = False
    dps1_burn_done = False
    burn1_done = False
    
    while current_time < total_duration:
        #Task 3: Include a retrograde rocket burn at 101104 seconds that reduces the velocity by 7.04 m/s
        # and include a rocket burn that increases the velocity at 212100 seconds by the amount given in the variable called boost.

        #Task 4: Implement Heun's method with adaptive step size. Note that the time is advanced at the end of this while loop.
        ###Your code here.
        acc = acceleration(current_time,position);
        xe = position + h*velocity
        ve = velocity + h*acc
        acce = acceleration(current_time+h, xe)
        
        position = position + h*0.5*(ve+velocity)
        velocity = velocity + h*0.5*(acc+acce)
        
        lte = numpy.linalg.norm(xe - position) + total_duration*numpy.linalg.norm(ve - velocity)
        h_new = h * math.sqrt( tolerance/lte )
        ###Your code here.

        h_new = min(0.5 * marker_time, max(0.1, h_new)) # restrict step size to reasonable range
            
        current_time += h
        h = h_new
        print current_time
        position_list.append(position.copy())
        velocity_list.append(velocity.copy())
        times_list.append(current_time)

    return position_list, velocity_list, times_list, boost

position, velocity, current_time, boost = apply_boost()

@show_plot
def plot_path(position_list, times_list):
    axes = matplotlib.pyplot.gca()
    axes.set_xlabel('Longitudinal position in m')
    axes.set_ylabel('Lateral position in m')
    previous_marker_number = -1;
    for position, current_time in zip(position_list, times_list):
         if current_time >= marker_time * previous_marker_number:
            previous_marker_number += 1
            matplotlib.pyplot.scatter(position[0], position[1], s = 2., facecolor = 'r', edgecolor = 'none')
            moon_pos = moon_position(current_time)
            if numpy.linalg.norm(position - moon_pos) < 30. * moon_radius: 
                axes.add_line(matplotlib.lines.Line2D([position[0], moon_pos[0]], [position[1], moon_pos[1]], alpha = 0.3, c = 'g')) 
    axes.add_patch(matplotlib.patches.CirclePolygon((0., 0.), earth_radius, facecolor = 'none', edgecolor = 'b'))
    for i in range(int(total_duration / marker_time)):
        moon_pos = moon_position(i * marker_time)
        axes.add_patch(matplotlib.patches.CirclePolygon(moon_pos, moon_radius, facecolor = 'none', edgecolor = 'g', alpha = 0.7))

    matplotlib.pyplot.axis('equal')

plot_path(position, current_time)