Source code for astroFunctions

#
#  ISC License
#
#  Copyright (c) 2016, Autonomous Vehicle Systems Lab, University of Colorado at Boulder
#
#  Permission to use, copy, modify, and/or distribute this software for any
#  purpose with or without fee is hereby granted, provided that the above
#  copyright notice and this permission notice appear in all copies.
#
#  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
#  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
#  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
#  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
#  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
#  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
#  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#


import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as la
from Basilisk.utilities import deprecated
from Basilisk.architecture import astroConstants

# Reference solar flux at Earth [W/m^2]
solarFluxEarth = 1372.5398

# ------------------------------------------------------------------------------------------------------------------- #
#                                           Normalize vector
#


[docs] def normalize(v): """normalize a vector""" norm = np.linalg.norm(v) if norm == 0: return v return v/norm
# ------------------------------------------------------------------------------------------------------------------- # # Orbital Period #
[docs] def orbitalPeriod(a, mu): """Return the orbit period""" P = 2*astroConstants.MPI * np.sqrt(a*a*a / mu) return P
[docs] def orbitalSMA(P, mu): """Return the semi-major axis""" a3 = mu * np.power(P / (2 * astroConstants.MPI), 2) exp = 1./3 a = np.power(a3, exp) return a
# ------------------------------------------------------------------------------------------------------------------- # # FROM OE 2 RV #
[docs] def Earth_RV(JDE): """return Earth (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'EARTH') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Mars_RV(JDE): """return Mars (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'MARS') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Jupiter_RV(JDE): """return Jupiter (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'JUPITER') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Venus_RV(JDE): """return Venus (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'VENUS') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Pluto_RV(JDE): """return Pluto (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'PLUTO') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Uranus_RV(JDE): """return Uranus (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'URANUS') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Neptune_RV(JDE): """return Neptune (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'NEPTUNE') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def Saturn_RV(JDE): """return Saturn (r,v)""" (a, e, i, Omega, w, nu) = ephemeridesMeeus(JDE, 'SATURN') return OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu)
[docs] def OE2RV(mu, a, e, i, Omega, w, nu): """OE to (r,v) conversion""" if e != 1: p = a*(1-e*e) else: print('ERROR: parabolic case') return c = np.cos(nu) s = np.sin(nu) r_PQW = np.array([p*c / (1 + e*c), p*s / (1 + e*c), 0]) v_PQW = np.array([-s * np.sqrt(mu/p), (e + c)*np.sqrt(mu/p), 0]) return PQW2IJK(Omega, i, w, r_PQW, v_PQW)
def quadrant4(x): return 2*astroConstants.MPI - x
[docs] def RV2OE(mu, r_IJK, v_IJK): """(r,v) to OE conversion""" r = r_IJK v = v_IJK K = np.array([0, 0, 1]) h = np.cross(r, v) n = np.cross(K, h) hm = np.linalg.norm(h) nm = np.linalg.norm(n) vm = np.linalg.norm(v) rm = np.linalg.norm(r) e_vec = 1./mu * ((vm*vm - mu/rm)*r - np.dot(r, v)*v) e = np.linalg.norm(e_vec) eta = 0.5*vm*vm - mu/rm if (e != 1.0): a = -0.5*mu/eta p = a*(1-e*e) else: a = 0. p = h*h/mu if h[2] == 0: i = 0.0 else: i = np.arccos(h[2]/hm) if n[0] == 0: Omega = 0.0 else: Omega = np.arccos(n[0]/nm) if n[1] < 0: Omega = quadrant4(Omega) dotProd = np.dot(n, e_vec) if dotProd == 0: omega = 0.0 else: omega = np.arccos(np.dot(n, e_vec) / (nm * e)) if e_vec[2] < 0: omega = quadrant4(omega) dotProd = np.dot(e_vec, r) if dotProd == 0: nu = 0.0 else: nu = np.arccos(np.dot(e_vec, r)/(e*rm)) if np.dot(r, v) < 0: nu = quadrant4(nu) return (a, e, i, Omega, omega, nu)
def PQW2IJK(Omega, i, w, r_PQW, v_PQW): sO = np.sin(Omega) cO = np.cos(Omega) si = np.sin(i) ci = np.cos(i) sw = np.sin(w) cw = np.cos(w) C11 = cO*cw - sO*sw*ci C12 = -cO*sw - sO*cw*ci C13 = sO*si C21 = sO*cw + cO*sw*ci C22 = -sO*sw + cO*cw*ci C23 = -cO*si C31 = sw*si C32 = cw*si C33 = ci C = np.array([[C11, C12, C13], [C21, C22, C23], [C31, C32, C33]]) r_IJK = C.dot(r_PQW) v_IJK = C.dot(v_PQW) return (r_IJK, v_IJK) # ------------------------------------------------------------------------------------------------------------------- # # MEEUS ALGORITHM # def ephemeridesMeeus(JDE, celestialBody): T = (JDE - 2451545.0)/36525 def ephem(vec, T): val = vec[0] + vec[1]*T + vec[2]*T + vec[3]*T*T*T return val if celestialBody == 'EARTH': L = 100.466449 + 35999.3728519 * T + (-0.00000568)*T*T + 0*T*T*T a = 1.000001018 e = 0.01670862 + (-0.000042037)*T + (-0.0000001236) * \ T*T + 0.00000000004*T*T*T i = 0 + 0.0130546*T + (-0.00000931)*T*T + (-0.000000034)*T*T*T Omega = 174.873174 + (-0.2410908)*T + 0.00004067 * \ T*T + (-0.000001327)*T*T*T Pi = 102.937348 + 0.3225557*T + 0.00015026*T*T + 0.000000478*T*T*T elif celestialBody == 'MARS': L = 355.433275 + 19140.2993313 * T + \ 0.00000261*T*T + (-0.000000003)*T*T*T a = 1.523679342 e = 0.09340062 + 0.000090483 * T + \ (-0.0000000806)*T*T + (-0.00000000035)*T*T*T i = 1.849726 + (-0.0081479) * T + (-0.00002255) * \ T*T + (-0.000000027)*T*T*T Omega = 49.558093 + (-0.2949846) * T + \ (-0.00063993)*T*T + (-0.000002143)*T*T*T Pi = 336.060234 + 0.4438898 * T + \ (- 0.00017321)*T*T + 0.000000300*T*T*T elif celestialBody == 'JUPITER': L = 34.351484 + 3034.9056746*T + (-0.00008501)*T*T + 0.000000004*T*T*T a = 5.202603191 + 0.0000001913*T e = 0.04849485 + 0.000163244*T + - \ 0.0000004719*T*T + (-0.00000000197)*T*T*T i = 1.303270 + (-0.0019872)*T + 0.00003318*T*T + 0.000000092*T*T*T Omega = 100.464441 + 0.1766828*T + \ 0.00090387*T*T + (-0.000007032)*T*T*T Pi = 14.331309 + 0.2155525*T + 0.00072252*T*T + (-0.000004590)*T*T*T elif celestialBody == 'VENUS': L = 181.979801 + 58517.8156760*T + \ 0.00000165*T*T + (-0.000000002)*T*T*T a = 0.72332982 e = 0.00677188 + (-0.000047766)*T + 0.0000000975 * \ T*T + 0.00000000044*T*T*T i = 3.394662 + (-0.0008568)*T + (-0.00003244)*T*T + 0.000000010*T*T*T Omega = 76.679920 + (-0.2780080)*T + (-0.00014256) * \ T*T + (-0.000000198)*T*T*T Pi = 131.563707 + 0.0048646*T + \ (-0.00138232)*T*T + (-0.000005332)*T*T*T elif celestialBody == 'PLUTO': L = 238.92903833 + 145.20780515*T a = 39.48211675 + (-0.00031596)*T e = 0.24882730 + 0.00005170*T i = 17.14001206 + 0.00004818*T Omega = 110.30393684 + (-0.01183482)*T Pi = 224.06891629 + (-0.04062942)*T elif celestialBody == 'SATURN': L_vec = np.array([50.077471, 1222.1137943, 0.00021004, -0.000000019]) a_vec = np.array([9.554909596, -0.0000021389, 0.0, 0.0]) e_vec = np.array([0.05550862, -0.000346818, - 0.0000006456, 0.00000000338]) i_vec = np.array([2.488878, 0.0025515, -0.00004903, 0.000000018]) Omega_vec = np.array( [113.665524, -0.2566649, -0.00018345, 0.000000357]) Pi_vec = np.array([93.056787, 0.5665496, 0.00052809, 0.000004882]) L = ephem(L_vec, T) a = ephem(a_vec, T) e = ephem(e_vec, T) i = ephem(i_vec, T) Omega = ephem(Omega_vec, T) Pi = ephem(Pi_vec, T) elif celestialBody == 'URANUS': L_vec = np.array([314.055005, 429.8640561, 0.00030434, 0.000000026]) a_vec = np.array([19.218446062, -0.0000000372, 0.00000000098, 0.0]) e_vec = np.array( [0.04629590, -0.000027337, 0.0000000790, 0.00000000025]) i_vec = np.array([0.773196, 0.0007744, 0.00003749, -0.000000092]) Omega_vec = np.array([74.005947, 0.5211258, 0.00133982, 0.000018516]) Pi_vec = np.array([173.005159, 1.4863784, 0.0021450, 0.000000433]) L = ephem(L_vec, T) a = ephem(a_vec, T) e = ephem(e_vec, T) i = ephem(i_vec, T) Omega = ephem(Omega_vec, T) Pi = ephem(Pi_vec, T) elif celestialBody == 'NEPTUNE': L_vec = np.array([304.348665, 219.8833092, 0.00030926, 0.000000018]) a_vec = np.array([30.110386869, -0.0000001663, 0.00000000069, 0.0]) e_vec = np.array([0.00898809, 0.000006408, - 0.0000000008, -0.00000000005]) i_vec = np.array([1.769952, -0.0093082, -0.00000708, 0.000000028]) Omega_vec = np.array([131.784057, 1.1022057, 0.00026006, -0.000000636]) Pi_vec = np.array([48.123691, 1.4262677, 0.00037918, -0.000000003]) L = ephem(L_vec, T) a = ephem(a_vec, T) e = ephem(e_vec, T) i = ephem(i_vec, T) Omega = ephem(Omega_vec, T) Pi = ephem(Pi_vec, T) else: print("Meeus coefficients for " + celestialBody + " not defined") L = 0. a = 0. e = 0. i = 0. Omega = 0. Pi = 0. return computeCOE(L, a, e, i, Omega, Pi) def computeCOE(L, a, e, i, Omega, Pi): # Units: L = L*astroConstants.D2R a = a*astroConstants.AU i = i*astroConstants.D2R Omega = Omega*astroConstants.D2R Pi = Pi*astroConstants.D2R w = Pi - Omega M = L - Pi C_cen = (2*e - 1//4 * np.power(e, 3) + 5//96 * np.power(e, 5)) * np.sin(M) +\ (5//4 * np.power(e, 2) - 11//24 * np.power(e, 4)) * np.sin(2*M) +\ (13//12 * np.power(e, 3) - 43//64 * np.power(e, 5)) * np.sin(3*M) +\ 103//96 * np.power(e, 4) * np.sin(4*M) +\ 1097//960 * np.power(e, 5) * np.sin(5*M) nu = M + C_cen return (a, e, i, Omega, w, nu) # ------------------------------------------------------------------------------------------------------------------- # # DATES CONVERSION # From GD to JD def JulianDate(GDE): yr = GDE[0] mo = GDE[1] d = GDE[2] JD = 367*yr - int(7*(yr + int((mo+9)/12))/4) + \ int(275*mo/9) + d + 1721013.5 return np.ceil(JD) # From JD to GD def GregorianDate(JDE): LMonth = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) T = (JDE - 2415019.5) / 365.25 Yr = 1900 + np.trunc(T) LeapYrs = np.trunc((Yr - 1900 - 1)*0.25) Days = JDE - 2415019.5 - ((Yr-1900)*365.0 + LeapYrs) if Days < 1.0: Yr = Yr - 1 LeapYrs = np.trunc((Yr - 1900 - 1)*0.25) Days = JDE - 2415019.5 - ((Yr-1900)*365.0 + LeapYrs) if Yr % 4 == 0: LMonth[1] = 29 DayOfYr = np.trunc(Days) sum = 0 i = 0 for month in LMonth: if (sum + month <= DayOfYr): sum += month i += 1 Mo = i + 1 Day = DayOfYr - sum GD = np.array([int(Yr), int(Mo), int(Day)]) return GD def exactGregorianDate(JD): L_month = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) T_1900 = (JD - 2415019.5)/365.25 Year = 1900 + np.trunc(T_1900) LeapYrs = np.trunc(0.25*(Year - 1900 - 1)) Days = (JD - 2415019.5) - (365.0*(Year - 1900) + LeapYrs) if Days < 1.0: Year = Year - 1 LeapYrs = np.trunc(0.25*(Year - 1900 - 1)) Days = (JD - 2415019.5) - (365.0*(Year - 1900) + LeapYrs) if Year % 4 == 0: L_month[1] = 29 DOY = np.trunc(Days) # Day of the year monthDays = 0 # summation of days in months m = 0 # months counter while True: if monthDays + L_month[m] >= DOY: break else: monthDays += L_month[m] m += 1 Month = m Day = DOY - monthDays time = (Days - DOY)*24 h = np.trunc(time) min = np.trunc(60*(time - h)) s = 3600 * (time - h - min/60.0) GD = '[' + str(int(Year)) + ', ' + str(int(Month) + 1) + ', ' + str(int(Day)) + ', ' + \ str(int(h)) + ':' + str(int(min)) + ':' + str(int(s)) + ']' return GD def optimalDate(GD0, DaysPastDeparture, TOF): JD0 = JulianDate(GD0) + DaysPastDeparture print('Departure Date: ', GregorianDate(JD0)) print('Arrival Date: ', GregorianDate(JD0 + TOF)) return (JD0, JD0+TOF) # ------------------------------------------------------------------------------------------------------------------- # # B PARAMETERS # def B_params_1(r, v, mu): rm = np.linalg.norm(r) vm = np.linalg.norm(v) h = np.cross(r, v) h_hat = normalize(h) e_vec = 1. / mu * ((vm * vm - mu / rm) * r - np.dot(r, v) * v) e = np.linalg.norm(e_vec) Betta = np.arccos(1. / e) a = 1. / (2 / rm - vm * vm / mu) b = np.abs(a) * np.sqrt(e * e - 1) k_hat = np.array([0, 0, 1]) S_hat = np.cos(Betta)/e * e_vec + np.sin(Betta) * \ normalize(np.cross(h_hat, e_vec)) T_hat = normalize(np.cross(S_hat, k_hat)) R_hat = np.cross(S_hat, T_hat) B_hat = np.cross(S_hat, h_hat) B = b*B_hat BT = np.dot(B, T_hat) BR = np.dot(B, R_hat) return (BT, BR) def B_params_2(v_in, v_out, mu): vm_in = np.linalg.norm(v_in) vm_out = np.linalg.norm(v_out) S = normalize(v_in) h = normalize(np.cross(v_in, v_out)) B_hat = normalize(np.cross(S, h)) k = np.array([0, 0, 1]) T = normalize(np.cross(S, k)) R = normalize(np.cross(S, T)) psi = np.arccos(np.dot(v_in, v_out) / (vm_in*vm_out)) rp = mu / (vm_in * vm_in) * \ (1.0 / np.cos(0.5*(astroConstants.MPI - psi)) - 1) temp = 1 + vm_in * vm_in * rp / mu B = mu / (vm_in * vm_in) * np.sqrt(temp*temp - 1) BT = B * np.dot(B_hat, T) BR = B * np.dot(B_hat, R) theta = np.arctan2(BR, BT) return (rp, psi*astroConstants.R2D, B, theta*astroConstants.R2D, BT, BR) # ------------------------------------------------------------------------------------------------------------------- # # RESONANT ORBITS # # v_inf_in: sc hyperbolic velocity # V1: Earth velocity at arrival (EGA1) # R1: Earth position at arrival (EGA1) # V2: Earth velocity at departure (EGA2) # N: Number of resonances def EarthResonantOrbit(v_inf_in, V1, R1, v_inf_out, V2, N): # compute P, a, v_scSun P = N * 365.242189 * 86400 # Earth orbital period n = P / (2.0*astroConstants.MPI) exp = 1./3 # Earth heliocentric orbit SMA a = np.power(n * n * astroConstants.MU_SUN, exp) rm_1 = la.norm(R1) vm_sc_h = np.sqrt(astroConstants.MU_SUN * (2./rm_1 - 1./a) ) # s/c heliocentric velocity # compute theta vm_inf_in = la.norm(v_inf_in) vm_P = la.norm(V1) cosTheta = (vm_inf_in * vm_inf_in + vm_P * vm_P - vm_sc_h * vm_sc_h) / (2 * vm_inf_in * vm_P) theta = np.arccos(cosTheta) # Compute T_vnc: DCM from VNC to ecliptic V_hat = normalize(V1) N_hat = normalize(np.cross(R1, V1)) C_hat = np.cross(V_hat, N_hat) T_vnc = np.array([V_hat, N_hat, C_hat]).T # Arrays to stack successful resonances that don't impact w/ the planet rp_GA1 = np.array([]) PHI_GA1 = np.array([]) rp_GA2 = np.array([]) PHI_GA2 = np.array([]) # Loop vm_inf_out = la.norm(v_inf_out) c = np.cos(astroConstants.MPI - theta) s = np.sin(astroConstants.MPI - theta) p = 360 # number of points phi_vec = np.linspace(0., 2*astroConstants.MPI, p) for phi in phi_vec: vGA_out_VNC = vm_inf_out * np.array([c, s*np.cos(phi), -s*np.sin(phi)]) v_GA1_out_ecl = np.dot(T_vnc, vGA_out_VNC) vm_GA1_out = la.norm(v_GA1_out_ecl) (rp, psi, B, theta, BT, BR) = B_params_2( v_inf_in, v_GA1_out_ecl, astroConstants.MU_EARTH) if rp > astroConstants.REQ_EARTH or rp < astroConstants.REQ_EARTH: rp_GA1 = np.append(rp_GA1, rp) PHI_GA1 = np.append(PHI_GA1, phi) v_GA2_in_ecl = v_GA1_out_ecl + V1 - V2 vm_GA2_in = la.norm(v_GA2_in_ecl) (rp, psi, B, theta, BT, BR) = B_params_2( v_GA2_in_ecl, v_inf_out, astroConstants.MU_EARTH) if rp > astroConstants.REQ_EARTH or rp < astroConstants.REQ_EARTH: rp_GA2 = np.append(rp_GA2, rp) PHI_GA2 = np.append(PHI_GA2, phi) plt.plot(PHI_GA1, rp_GA1, 'm', PHI_GA2, rp_GA2, 'b') plt.axhline(astroConstants.REQ_EARTH, color='k') plt.legend(['EGA 1', 'EGA 2', 'r$_{P, min}$'], loc='lower right') plt.xlabel(r'$\phi$ [rad]') plt.ylabel('Perigee Radius [km]') # plt.show() return def testResonance(): R1 = np.array([31161153575.0634, 143995536213.184, 3018421.6823707]) R2 = np.array([31194146713.798, 143988869052.709, 2352102.45497072]) V1 = np.array([ -29599.9908905588, 6189.32095554111, 0.0729005972949417]) V2 = np.array([-29598.5134550237, 6195.92283490717, 0.0568758149973483]) v_inf_in = np.array([-0.320955575422136, -8.79065711515436, 1.70894261814904]) v_inf_out = np.array([-8.95399388012787, -0.346120703226863, -0.350867778719334]) N = 2 EarthResonantOrbit(v_inf_in, V1/1000, R1/1000, v_inf_out, V2/1000, N) # ------------------------------------------------------------------------------------------------------------------- # # HELIOCENTRIC ORBIT PLOT # def plotPlanetOrbit(JD, planet, color, ax): rx = np.array([]) ry = np.array([]) rz = np.array([]) (a, e, i, Omega, w, nu0) = ephemeridesMeeus(JD, planet) nu_vec = np.linspace(nu0, nu0 + 2*astroConstants.MPI, 360) for nu in nu_vec: (r, v) = OE2RV(astroConstants.MU_SUN, a, e, i, Omega, w, nu) rx = np.append(rx, r[0]) ry = np.append(ry, r[1]) rz = np.append(rz, r[2]) ax.plot(rx, ry, rz, c=color) ax.scatter(rx[0], ry[0], rz[0], c=color) max_range = np.array([rx.max()-rx.min(), ry.max() - ry.min(), rz.max()-rz.min()]).max() Xb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, - 1:2:2][0].flatten() + 0.5*(rx.max()+rx.min()) Yb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, - 1:2:2][1].flatten() + 0.5*(ry.max()+ry.min()) Zb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, - 1:2:2][2].flatten() + 0.5*(rz.max()+rz.min()) # Comment or uncomment following both lines to test the fake bounding box: for xb, yb, zb in zip(Xb, Yb, Zb): ax.plot([xb], [yb], [zb], 'w') def plotSolarSystem(JD, JD_ref): fig = plt.figure() ax = fig.add_subplot(111, projection='3d') plotPlanetOrbit(JD, 'EARTH', 'dodgerblue', ax) plotPlanetOrbit(JD, 'SATURN', 'orchid', ax) plotPlanetOrbit(JD, 'URANUS', 'lightgreen', ax) plt.title('Days Past JD_ref = ' + str(JD - JD_ref)) plt.show() # ------------------------------------------------------------------------------------------------------------------- # # TISSERAND PLOT # def V_circular(r, mu): return np.sqrt(mu / r) def rotationMatrix(alpha): M = np.identity(3) M[0, 0] = np.cos(alpha) M[0, 1] = -np.sin(alpha) M[1, 0] = np.sin(alpha) M[0, 1] = np.cos(alpha) return M def TisserandPlot(R_P, v_inf, color): RP = np.array([]) RA = np.array([]) R_P_vec = np.array([0., - R_P, 0.]) V_P = V_circular(R_P, astroConstants.MU_SUN) V_P_vec = np.array([V_P, 0., 0.]) points = 360 alpha_vec = np.linspace(0., astroConstants.MPI, points) for alpha in alpha_vec: M = rotationMatrix(alpha) V_P_hat = normalize(V_P_vec) v_inf_vec = v_inf * np.dot(M, V_P_hat) V_inf_vec = V_P_vec + v_inf_vec (a, e, i, Omega, omega, nu) = RV2OE( astroConstants.MU_SUN, R_P_vec, V_inf_vec) rp = a * (1. - e) ra = a * (1. + e) RP = np.append(RP, rp/astroConstants.AU) RA = np.append(RA, ra/astroConstants.AU) plt.loglog(RP, RA, color) def Tisserand_P1_2_P2(r_planet_vec, v_inf_P1_vec, v_inf_P2_vec): for r_planet in r_planet_vec: plt.loglog(r_planet/astroConstants.AU, r_planet/astroConstants.AU, 'ko') color_vec = np.array( ['dodgerblue', 'lightgreen', 'magenta', 'lightsalmon', 'indigo']) color_vec = np.array(['k', 'k', 'k', 'k', 'k']) i = 0 for v_inf in v_inf_P1_vec: TisserandPlot(r_planet_vec[0], v_inf, color_vec[i]) i += 1 j = 0 for v_inf in v_inf_P2_vec: TisserandPlot(r_planet_vec[1], v_inf, color_vec[j]) j += 1 # x-axis plot limits l = len(r_planet_vec) - 1 rp_min = r_planet_vec[0] / astroConstants.AU * 0.2 rp_max = r_planet_vec[l] / astroConstants.AU * 1.8 # y-axis plot limits ra_min = r_planet_vec[0] / astroConstants.AU * 0.2 ra_max = r_planet_vec[l] / astroConstants.AU * 1.8 plt.xlim([rp_min, rp_max]) plt.ylim([ra_min, ra_max]) plt.xlabel('Radius of Periapse [AU]') plt.ylabel('Radius of Apoapse [AU]') plt.title(r'$v_{\infty}^{Earth} $ = ' + str(v_inf_P1_vec) + r' km/s \n' r'$v_{\infty}^{Saturn} $ = ' + str(v_inf_P2_vec) + ' km/s') def Tisserand_ESU(): r_planet_vec = np.array( [astroConstants.SMA_EARTH, astroConstants.SMA_SATURN, astroConstants.SMA_URANUS]) v_inf_E_vec = np.array([9, 11, 13, 15, 17]) v_inf_S_vec = np.array([6, 7, 8, 9, 10]) Tisserand_P1_2_P2(r_planet_vec, v_inf_E_vec, v_inf_S_vec) TisserandPlot(astroConstants.SMA_EARTH, v_inf_E_vec[2], 'cyan') TisserandPlot(astroConstants.SMA_SATURN, v_inf_S_vec[3], 'cyan') plt.show() # ------------------------------------------------------------------------------------------------------------------- # # DEPRECATED CONSTANTS # def __getattr__(name): """ Special `__getattr__` to handle deprecation warning """ AU = 1.49597870691e8 # [km] astroConstants = { # Useful conversions 'R2D': 180.0 / np.pi, # [deg] 'D2R': np.pi / 180.0, # [rad] 'AU': AU, 'SEC2DAY': 1.0 / (3600*24), 'DAY2SEC': 1.0 * 3600 * 24, 'PI': np.pi, # Gravitational parameters [km^3 / s^2] 'mu_S': 1.32712440018e11, 'mu_E': 3.986004415e5, 'mu_M': 4.28283100e4, 'mu_V': 3.24858599e5, 'mu_J': 1.26686534e8, 'mu_Saturn': 3.79395000e7, 'mu_U': 5.79396566e6, 'mu_N': 6.83509920e6, 'mu_Moon': 4902.799, # Planets' Radius [km] 'S_radius': 695508, 'J_radius': 71492, 'V_radius': 6051.8, 'E_radius': 6378.1363, 'M_radius': 3396.19, 'Saturn_radius': 60268, 'U_radius': 25559, 'N_radius': 24764, 'Moon_radius': 1738.1, # Planetary Heliocentric orbits [km] 'a_E': 1.0 * AU, 'a_M': 1.52367934 * AU, 'a_Mercury': 0.387 * AU, 'a_V': 0.723 * AU, 'a_J': 5.203 * AU, 'a_Saturn': 9.537 * AU, 'a_U': 19.191 * AU, 'a_N': 30.069 * AU, 'a_P': 39.482 * AU, } if name in astroConstants: deprecated.deprecationWarn( name, '2025/08/31', 'Astro constants should be used from Basilisk.architecture.astroConstants') return astroConstants[name] def main(): Mars_RV(JulianDate([2018, 10, 16])) if __name__ == '__main__': main()