ballistic_sample2

Using ballistic module to simulate Moon orbiter

""" Using ballistic module to simulate Moon orbiter
"""
import numpy as np
import matplotlib.pyplot as plt

from irfpy.util import ballistic
from irfpy.util.planets import moon

def main():

    msc = 700.0      # It does not matter, but anyway spacecraft mass
    m = moon['mass']
    rm = 1738e3

    ### Spacecraft orbit determination
    d_apo = 270 * 1e3 + rm    # 270 km above the surface
    d_per = 30 * 1e3 + rm    # 30 km above the surface

    v_per, v_apo = ballistic.speeds_from_distances(d_per, d_apo, m)

    print('APO: d={:.3f} km : v={:.3f} km/s : w={:.3e} rad/s'.format(d_apo, v_apo, v_apo / d_apo))
    print('PER: d={:.3f} km : v={:.3f} km/s : w={:.3e} rad/s'.format(d_per, v_per, v_per / d_per))

    ### Spacecraft initial state (apocenter at north pole)
    r0 = d_apo
    th0 = np.deg2rad(70)    # North pole.
    v0 = 0       # Pure horizontal
    w0 = v_apo / d_apo     # Omega

    ### Ballistic trajectory
    traj = ballistic.Trajectory(r0, th0, v0, w0, m, m=msc)
#    traj.print_values()

    ### theta dependence
    theta_list = np.linspace(0, 360, 361)
    tlist = np.deg2rad(theta_list)
    rlist = []
    vlist = []
    wlist = []

    for theta in tlist:
        r = traj.r(theta)
        v = traj.v(theta)
        w = traj.omega(theta)

        rlist.append(r)
        vlist.append(v)
        wlist.append(w)

    rlist = np.array(rlist)
    vlist = np.array(vlist)
    wlist = np.array(wlist)

    ### Plotting

    plt.figure()
    plt.subplot(4, 1, 1)
    plt.plot(theta_list, (rlist - rm)/1000, 'b')
#    plt.axhline(rm, c='y')
    plt.ylabel('h [km]')
    plt.subplot(4, 1, 2)
    plt.plot(theta_list, vlist, 'b')
    plt.axhline(0, c='k')
    plt.ylabel('vr [m/s]')
    plt.subplot(4, 1, 3)
    plt.plot(theta_list, wlist * rlist, 'b')
    plt.axhline(0, c='k')
    plt.ylabel('vt [m/s]')
    plt.subplot(4, 1, 4)
    plt.plot(theta_list, np.rad2deg(np.arctan2(vlist, wlist * rlist)), 'b')
    plt.axhline(0, c='k')
    plt.axvline(270, c='k', linestyle=":")
    plt.ylabel('v_angle [deg]')

    plt.figure()
    plt.subplot(111)

    theta_list_moon = np.linspace(0, 360, 361)
    tlistmoon = np.deg2rad(theta_list_moon)
    xmoon = rm * np.cos(tlistmoon)
    ymoon = rm * np.sin(tlistmoon)

    xlist = rlist * np.cos(tlist)
    ylist = rlist * np.sin(tlist)

    plt.plot(xlist, ylist, 'b')
    plt.plot(xmoon, ymoon, 'y')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect(1)

    plt.show()

if __name__ == '__main__':
    main()