ballistic_sample

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(90)    # 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(90, 270, 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)

    ### Impact probe speed at the south pole, first delta-v.
    second_dv_angle = 540

    thip0 = 3 * np.pi / 2.
    rip0 = traj.r(thip0)
    vip0 = traj.v(thip0)
    wip0 = traj.omega(thip0)

    delta_vt = 30
    delta_vr = 0

    vip0 += delta_vr
    vtip0 = rip0 * wip0    # Speed
    vtip0 += delta_vt    # deltav
    wip0 = vtip0 / rip0   # Back to angular speed

    mip = 10   # kg

    trajip = ballistic.Trajectory(rip0, thip0, vip0, wip0, m, m=mip)

    theta_listip = np.linspace(270, second_dv_angle, 181)
    tlistip = np.deg2rad(theta_listip)
    riplist = []
    viplist = []
    wiplist = []

    for theta in tlistip:
        r = trajip.r(theta)
        v = trajip.v(theta)
        w = trajip.omega(theta)

        riplist.append(r)
        viplist.append(v)
        wiplist.append(w)

    riplist = np.array(riplist)
    viplist = np.array(viplist)
    wiplist = np.array(wiplist)

    ### 2nd delta-v
    thip1 = np.deg2rad(second_dv_angle)
    rip1 = riplist[-1]
    vip1 = viplist[-1]
    wip1 = wiplist[-1]

    delta_vr1 = 0
    delta_vt1 = -10
    vip1 += delta_vr1
    vtip1 = rip1 * wip1    # Speed
    vtip1 += delta_vt1    # deltav
    wip1 = vtip1 / rip1   # Back to angular speed


    trajip1 = ballistic.Trajectory(rip1, thip1, vip1, wip1, m, m=mip)

    theta_listip1 = np.linspace(second_dv_angle, 720, 1801)
    tlistip1 = np.deg2rad(theta_listip1)
    riplist1 = []
    viplist1 = []
    wiplist1 = []

    for theta in tlistip1:
        r = trajip1.r(theta)
        v = trajip1.v(theta)
        w = trajip1.omega(theta)

        riplist1.append(r)
        viplist1.append(v)
        wiplist1.append(w)

    riplist1 = np.array(riplist1)
    viplist1 = np.array(viplist1)
    wiplist1 = np.array(wiplist1)


    ### Plotting

    plt.figure()
    plt.subplot(4, 1, 1)
    plt.plot(theta_list, rlist, 'b')
    plt.plot(theta_listip, riplist, 'g')
    plt.plot(theta_listip1, riplist1, 'r')
    plt.axhline(rm, c='y')
    plt.ylabel('r [m]')
    plt.subplot(4, 1, 2)
    plt.plot(theta_list, vlist, 'b')
    plt.plot(theta_listip, viplist, 'g')
    plt.plot(theta_listip1, viplist1, 'r')
    plt.axhline(0, c='k')
    plt.ylabel('vr [m/s]')
    plt.subplot(4, 1, 3)
#    plt.plot(theta_list, wlist, 'b')
#    plt.plot(theta_listip, wiplist, 'g')
#    plt.plot(theta_listip1, wiplist1, 'r')
#    plt.ylabel('omega [rad/s]')
    plt.plot(theta_list, wlist * rlist, 'b')
    plt.plot(theta_listip, wiplist * riplist, 'g')
    plt.plot(theta_listip1, wiplist1 * riplist1, 'r')
    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.plot(theta_listip, np.rad2deg(np.arctan2(viplist, wiplist * riplist)), 'g')
    plt.plot(theta_listip1, np.rad2deg(np.arctan2(viplist1, wiplist1 * riplist1)), 'g')
    plt.axhline(0, c='k')
    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)

    xiplist = riplist * np.cos(tlistip)
    yiplist = riplist * np.sin(tlistip)

    xiplist1 = riplist1 * np.cos(tlistip1)
    yiplist1 = riplist1 * np.sin(tlistip1)

    plt.plot(xlist, ylist, 'b')
    plt.plot(xiplist, yiplist, 'g')
    plt.plot(xiplist1, yiplist1, 'r')
    plt.plot(xmoon, ymoon, 'y')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect(1)


    plt.show()

if __name__ == '__main__':
    main()