Перелет по Гоману

Рейтинг: 0Ответов: 1Опубликовано: 23.05.2023

У меня есть функция 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

Ответы

▲ 0

runtime система говорит 'KeyError'. И указывает оператор, на котором выявлено недопустимое значение ключа:

posCent, velCent = kernel[CENTR,cent].compute_and_differentiate(jd)

Вот и проверяйте значения переменных CENTR и cent - почему они вылазят за границы индексов массива...