# 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)
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)
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
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)
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.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)
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()