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