Перелет по Гоману
У меня есть функция Ephemeris
def Ephemeris(jd, targ, cent, UnitR, UnitV):
##_____EPHEM__________________________
#kernel = SPK.open('de405.bsp')
## print(kernel)
#posMoon, velMoon = kernel[3,301].compute_and_differentiate(2459580.5)#todaydate
#posEarth, velEarth = kernel[3,399].compute_and_differentiate(2459580.5)#todaydate
#position = numpy.zeros(3)
#velocity = numpy.zeros(3)
#for i in range(3):
# position[i]= posMoon[i]-posEarth[i]
# velocity[i]=(velMoon[i]-velEarth[i])/86400
#print("moon position = ", position, "km")
#print("moon velosity = ", velocity, "km/sec")
#Solar System Barycenter (0) -> Mercury Barycenter (1)
#Solar System Barycenter (0) -> Venus Barycenter (2)
#Solar System Barycenter (0) -> Earth Barycenter (3)
#Solar System Barycenter (0) -> Mars Barycenter (4)
#Solar System Barycenter (0) -> Jupiter Barycenter (5)
#Solar System Barycenter (0) -> Saturn Barycenter (6)
#Solar System Barycenter (0) -> Uranus Barycenter (7)
#Solar System Barycenter (0) -> Neptune Barycenter (8)
#Solar System Barycenter (0) -> Pluto Barycenter (9)
#Solar System Barycenter (0) -> Sun (10)
# Earth Barycenter (3) -> Moon (301)
# Earth Barycenter (3) -> Earth (399)
#____________________________________
RV = np.zeros(6);
# kernel = SPK.open('de721.bsp');
kernel = SPK.open('de405.bsp');
if 0 <= targ <= 11 and cent==0 or 300 < targ < 400 and cent==3:
posTarg, velTarg = kernel[cent,targ].compute_and_differentiate(jd);
for i in range(3):
RV[i] = posTarg[i]/UnitR;
RV[3+i] = velTarg[i]/86400/UnitV;
return RV;
if targ==301 or targ==399:
CENTR = 3;
posTarg, velTarg = kernel[CENTR,targ].compute_and_differentiate(jd);
elif 0 <= targ <= 11:
CENTR = 0;
posTarg, velTarg = kernel[CENTR,targ].compute_and_differentiate(jd);
if cent==301 or cent==399:
CENTR = 3;
posCent, velCent = kernel[CENTR,cent].compute_and_differentiate(jd);
elif 0 <= cent <= 11:
CENTR = 0;
posCent, velCent = kernel[CENTR,cent].compute_and_differentiate(jd);
if 0 <= targ <= 11 and 300 < cent < 400 or 300 < targ < 400 and 0 <= cent <= 11:
posCENTR, velCENTR = kernel[0,3].compute_and_differentiate(jd);
if 0 <= targ <= 11:
for i in range(3):
RV[i] = (posTarg[i]-posCent[i]-posCENTR[i])/UnitR;
RV[3+i] = (velTarg[i]-velCent[i]-velCENTR[i])/86400/UnitV;
elif 0 <= cent <= 11:
for i in range(3):
RV[i] = (posTarg[i]+posCENTR[i]-posCent[i])/UnitR;
RV[3+i] = (velTarg[i]+velCENTR[i]-velCent[i])/86400/UnitV;
else:
for i in range(3):
RV[i] = (posTarg[i]-posCent[i])/UnitR;
RV[3+i] = (velTarg[i]-velCent[i])/86400/UnitV;
return RV
И есть функция вычисляющая параметры гомана
def paramhelio(r1, r2, mu, year, month, day, Ep, Jp,phi_Jup,phi_E):
dv1 = np.sqrt(mu / r1) * (np.sqrt(2 * r2 / (r1 + r2)) - 1)
dv2 = np.sqrt(mu / r2) * (1 - np.sqrt(2 * r1 / (r1 + r2)))
v_total = abs(dv1) + abs(dv2)
flight_time = np.pi * np.sqrt((r1 + r2) ** 3 / (8 * mu))
dphih = np.pi*(1-np.sqrt(((r1/r2)+1)**3)/np.sqrt(8))
n1 = np.sqrt(1/1**3)
n2 = np.sqrt(1/5**3)
startdate = toJD(year,month,day) + ((np.pi - n2 * flight_time/31536000 - phi_Jup + phi_E) / abs(n2 - n1)) + (2 * np.pi / abs(n2-n1)) * i
return dv1, dv2, v_total, flight_time, dphih, toCAL(startdate)
У меня получилось нарисовать траекторию перелета от Земли до Юпитера
def plotH(r1, r2, year, month, day):
param = paramhelio(r1, r2, fM_Sun, year, month, day, Ep, Jp,phi_Jup,phi_E)
jd0 = toJD(year,month,day)
jdk = toJD(param[5][0],param[5][1],param[5][2])
EphE = Ephemeris(jd0, 3, 0, 1, 1)
EphJ = Ephemeris(jdk, 5, 0, 1, 1)
SE = np.sqrt(EphE[0]**2 + EphE[1]**2)
SJ = np.sqrt(EphJ[0]**2 + EphJ[1]**2)
R = 20000000
figure, axes = plt.subplots()
plt.axis([-1000000000, 1000000000, -1000000000, 1000000000])
axes.set_aspect(1)
c1 = plt.Circle((0,0), radius = R, color = "y")
c2 = plt.Circle((0,0), radius = SJ, fill = False)
c3 = plt.Circle((0,0), radius = SE, fill = False)
c4 = plt.Circle((EphE[0],EphE[1]), radius = R, color = "b")
c5 = plt.Circle((EphJ[0],EphJ[1]), radius = R, color = "g")
b = np.sqrt(r1*r2)
t = np.linspace(-90, 90, 180)
x = b * np.cos(np.radians(t))
y = (SE + SJ) / 2 * np.sin(np.radians(t)) + (SE + SJ) / 2
x_n = x * np.cos(arctan2(EphJ[1],EphJ[0]) - np.pi / 2) - y * np.sin(arctan2(EphJ[1],EphJ[0]) - np.pi / 2) + EphE[0]
y_n = y * np.cos(arctan2(EphJ[1],EphJ[0]) - np.pi / 2) - x * np.sin(arctan2(EphJ[1],EphJ[0]) - np.pi / 2) + EphE[1]
plt.plot(x_n,y_n, color="r")
plt.gca().add_artist(c1)
plt.gca().add_artist(c2)
plt.gca().add_artist(c3)
plt.gca().add_artist(c4)
plt.gca().add_artist(c5)
Теперь хочу нарисовать траекторию движения относительно Ио (спутник Юпитера), но когда хочу взять данные из Ephemeris для Ио:
EphIo = Ephemeris(jd, 301, 0, 1, 1)
мне выдает ошибку:
KeyError
Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_7572\1480923696.py in <module>
2 jd0 = toJD(year, month, day)
3 jdk = toJD(param[5][0], param[5][1], param[5][2])
----> 4 EphE = Ephemeris(jdk, 301, 0, 1, 1)
5 print(EphE)
~\AppData\Local\Temp\ipykernel_7572\3460855085.py in Ephemeris(jd, targ, cent, UnitR, UnitV)
54 elif 0 <= cent <= 11:
55 CENTR = 0;
---> 56 posCent, velCent = kernel[CENTR,cent].compute_and_differentiate(jd);
57
Источник: Stack Overflow на русском