Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import os
def omega_bar_dot(i, *data):
n, J2, R, a, e = data
return 3/4*n*J2*R**2/a**2*(5*np.cos(i)**2 - 1)/(1-e**2)**2
def main():
GM = 3.986*10**5 # km3/s2.
J2 = 0.00108
R = 6370 # km
rps = np.arange(600, 20307, 10) # km
T = 1/3*24*60*60 # s
# Get SMA
a = (GM * T ** 2 / (4 * np.pi ** 2)) ** (1 / 3) # km
es = 1 - rps / a
i = np.arccos(np.sqrt(1 / 5))
n = 2 * np.pi / T
Omega_bar_dots = -3 / 2 * n * J2 * (R / a) ** 2 * np.cos(i) / (1 - es ** 2) ** 2
# Find Omega_bar_dot where precession rate is 360 deg/year
Omega_bar_dot_req = 360*np.pi/180/(60*60*24*365)
# get index where where the omega_bar_dot are less than this
idx = np.where(abs(Omega_bar_dots) < Omega_bar_dot_req)[0][0]
# Choose Rp where this happens
rp_choose = rps[idx]
a = (GM * T ** 2 / (4 * np.pi ** 2)) ** (1 / 3) # km
e = 1 - rp_choose / a
i = np.arccos(np.sqrt(1 / 5))
i_neg = np.arccos(-np.sqrt(1 / 5))
n = 2 * np.pi / T
Omega_bar_dot = -3 / 2 * n * J2 * (R / a) ** 2 * np.cos(i) / (1 - e ** 2) ** 2
ra = a*(1+e)
print("Choose r_p [km]: ", rp_choose)
print("SMA [km]: ", a)
print("Eccentricity: ", e)
print("Inclination [deg]: ", np.rad2deg(i), " or ", np.rad2deg(i_neg))
print("RAAN rate [rad/s]: ", Omega_bar_dot, " or ", -Omega_bar_dot)
print("RAAN rate [deg/yr]: ", abs(Omega_bar_dot)*60*60*24*365*180/(np.pi))
print("r_a [km]: ", ra)
# Plot omega_bar_dot vs. i
i = np.arange(0, np.pi, 0.01)
omega_bar_dot_plot = omega_bar_dot(i, n, J2, R, a, e)
plt.plot(np.rad2deg(i), omega_bar_dot_plot)
plt.axhline(y=0, color='r')
plt.axvline(x = 63.43494882292201, color='r')
plt.axvline(x = 116.56505117707799, color='r')
plt.xlabel("Inclination [deg]")
plt.ylabel(r"$\dot{\bar{\Omega}}$ [rad/s]")
plt.grid()
plt.savefig(os.path.join("./figs", "incSols.png"))
# plot Rp and eccentricity vs. Omega_bar_dot
# Make subplot for Rp
fig, axs = plt.subplots(2)
fig.suptitle("Problem 1")
axs[0].plot(rps, Omega_bar_dots * 60 * 60 * 24 * 365 * 180 / (np.pi))
axs[0].set_xlabel("Radius of Periapsis [km]")
axs[0].set_ylabel(r"$\dot{\bar{\Omega}}$ [deg/yr]")
# make a vertical line at x = rp_choose
axs[0].axvline(x=rp_choose, color='r')
axs[0].grid()
axs[0].legend([r"$\dot{\bar{\Omega}}$", "r_p chosen"])
# make subplot for eccentricity versus Omega_bar_dot
axs[1].plot(es, Omega_bar_dots * 60 * 60 * 24 * 365 * 180 / (np.pi))
axs[1].set_xlabel("Eccentricity")
axs[1].set_ylabel(r"$\dot{\bar{\Omega}}$ [deg/yr]")
# make a vertical line at x = rp_choose
axs[1].axvline(x=e, color='r')
axs[1].grid()
axs[1].legend([r"$\dot{\bar{\Omega}}$", "e chosen"])
# Increase vertical padding between subplots
fig.tight_layout(pad=1.0)
plt.savefig(os.path.join("./figs", "prob1.png"))
plt.close()
if __name__ == "__main__":
main()