Попытка решить спутниковое уравнение, но я получаю пустой сюжет ⇐ Python

[Anonymous]

Я пытаюсь численно интегрировать уравнения движения для спутника, движущегося вокруг сфероидальной земли, используя метод Хеуна, но я получаю Runtimewarnning: недействительное значение, встречающееся в скалярном разделении, с последующим пустым графиком, где он показывает только оси без каких -либо кривых. Вот объяснение уравнений, которые я пытаюсь решить: < /p>

Я ожидал получить alliptical, как или изначально, что -то, что -то, что -то, похожее: что -то, что -то, что -то, что -то, похожее: что -то, что -то, похожее: что -то, что -то, что -то, что -то, похожее: что -то, что -то, похожее: что -то, что -то, похожее: что -то, что -то, что -то, похожее: что -то, что -то, что -то, похожее: что -то, что -то, что -то, похожее: что -то, что -то, что -то, похожее: что -то, что -то вроде: alt = "спутниковая орбита" src = "https://i.sstatic.net/inccj8j8.png"/>
Вот мой код:
from math import sin, cos, pi
import numpy as np
from matplotlib import pyplot as plt
G = 6.67e-11
M = 5.97e24
m = 3000
R = 6371000
J2 = 0.001341125
dt = 0.1
t = np.linspace(0,500,5000)
r = np.empty(len(t))
phi = np.empty(len(t))
th = np.empty(len(t))
p_r = np.empty(len(t))
p_phi = np.empty(len(t))
p_th = np.empty(len(t))

def f(r,phi,th,p_r,p_phi,p_th):
return p_r/m
def g(r,phi,th,p_r,p_phi,p_th):
return p_phi**2/(m*r**3) + p_th**2/(m*r**3*cos(phi)**2) - G*M*m/r**2 + 3*J2*G*M*m*R**2*(3*sin(phi)**2 - 1)/(2*r**4)
def h(r,phi,th,p_r,p_phi,p_th):
return p_phi/(m*r**2)
def j(r,phi,th,p_r,p_phi,p_th):
return -p_th**2*sin(phi)/(m*r**2*cos(phi)**3) - 3*J2*G*M*m*R**2*sin(phi)*cos(phi)/r**3
def k(r,phi,th,p_r,p_phi,p_th):
return p_th/(m*r**2*cos(phi)**2)
def Heun(f,g,h,j,k,t):
r[0] = 1.4371e7
phi[0] = 0.49672
th[0] = -1.4055
p_r[0] = 1.1775e7
p_phi[0] = 1.744e14
p_th[0] = 5.3921e14
rp = np.empty(len(t))
phip = np.empty(len(t))
thp = np.empty(len(t))
p_rp = np.empty(len(t))
p_phip = np.empty(len(t))
p_thp = np.empty(len(t))
for n in range(len(t) - 1):
p_rp[n+1] = p_r[n] + dt*g(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n])
rp[n+1] = r[n] + dt*p_r[n]/m
p_phip[n+1] = p_phi[n] + dt*j(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n])
phip[n+1] = phi[n] + dt*p_phi[n]/(m*r[n]**2)
p_thp[n+1] = p_th[n]
thp[n+1] = th[n] + dt*k(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n])
p_r[n+1] = p_r[n] + dt/2*(g(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n]) + g(rp[n],phip[n],thp[n],p_rp[n],p_phip[n],p_thp[n]))
r[n+1] = r[n] + dt/2*(p_r[n]/m + p_rp[n]/m)
p_phi[n+1] = p_phi[n] + dt/2*(j(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n]) + j(rp[n],phip[n],thp[n],p_rp[n],p_phip[n],p_thp[n]))
phi[n+1] = phi[n] + dt/2*(h(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n]) + h(rp[n],phip[n],thp[n],p_rp[n],p_phip[n],p_thp[n]))
p_th[n+1] = p_th[n]
th[n+1] = th[n] + dt/2*(k(r[n],phi[n],th[n],p_r[n],p_phi[n],p_th[n]) + k(rp[n],phip[n],thp[n],p_rp[n],p_phip[n],p_thp[n]))
if phi[n+1] == pi/2:
break
return t,r,phi,th,p_r,p_phi,p_th
t,r,phi,th,p_r,p_phi,p_th = Heun(f,g,h,j,k,t)
x = r*np.cos(phi)*np.cos(th)
y = r*np.cos(phi)*np.sin(th)
z = r*np.sin(phi)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(x, y, z, 'blue')
plt.show()



Подробнее здесь: https://stackoverflow.com/questions/797 ... empty-plot