轨道衰减仿真中实现表面碰撞时停止仿真的技术求助

我正在做一个轨道衰减仿真，遇到了个问题：当物体碰撞到行星表面时，仿真不会停止/中断，而是会继续在行星内部运行，没法精准获取碰撞发生的时刻。

我的代码如下： #

import numpy as np
import scipy.integrate as sci
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

plt.close()

G = 6.6742e-11

###Planets
###Earth
name = 'Earth'
Rplanet= 6357000 #m
mplanet= 5.972e24 #Kg

###Cube satellite
mass = 1 #Kg
S = 1/100 #Transversal area
CD = 2

###Initial Conditions
'''
I'll have a R3 cartesian system. Therefore i'll have to set the
planet in terms of xyz and Vector position, speed and force
in those unit vectors. So i'll define all vectors as arrays.
'''
x0 = 0
y0 = Rplanet + 350000
z0 = 0
r0 = Rplanet + 350000
velx0 = 0
vely0 = 0
velz0 = np.sqrt((G*mplanet)/r0)
period = 1e5

xdot = velx0
ydot = vely0
zdot = velz0

###AeroDynamics Model
class Aerodynamics():
    def __init__(self, name):
        self.name = name
        if name == 'Earth':
            ###Model with earth data
            self.beta = 0.1354/1000 #Density constant
            self.rhos = 1.225 #kg/m^3
    def getDensity(self, altitude):
        if self.name == 'Earth':
            rho = self.rhos*np.exp(-altitude*self.beta)
            return rho

aeroModel = Aerodynamics(name)



###Gravity model

def gravity(x, y, z):
    global Rplanet, mplanet, G

    r = np.sqrt(x**2 + y**2 + z**2)

    if r < Rplanet:
        accelx = 0
        accely = 0
        accelz = 0
    else:
        accelx = (G*mplanet)*x/(r**3)
        accely = (G*mplanet)*y/(r**3)
        accelz = (G*mplanet)*z/(r**3)
    return np.asarray([accelx, accely, accelz]), r
    


###Differential equation of the system
def Derivatives(state, t):
    global mass

    x = state[0]
    y = state[1]
    z = state[2]
    velx = state[3]
    vely = state[4]
    velz = state[5]

    xdot = velx
    ydot = vely
    zdot = velz

    ###Forces
    ###Gravity
    accel, r = gravity(x, y, z)
    GravityF = -accel*mass

    ###Aerodynamics
    altitude = r - Rplanet
    rho = aeroModel.getDensity(altitude)
    V = np.sqrt(velx**2 + vely**2 + velz**2)
    qinf = -1/2*rho*abs(V)*S
    aeroF = qinf*CD*np.asarray([velx, vely, velz])
    
    Forces = GravityF + aeroF

    ddot = Forces/mass

    statedot = np.asarray([xdot, ydot, zdot, ddot[0], ddot[1], ddot[2]])
    return statedot


###Run

stateInitial = np.asarray([x0, y0, z0, velx0, vely0, velz0])

###Time window

tout = np.linspace(0, period, 1000)

stateout = sci.odeint(Derivatives,stateInitial,tout)

xout = stateout[:,0]
yout = stateout[:,1]
zout = stateout[:,2]
velxout = stateout[:,3]
velyout = stateout[:,4]
velzout = stateout[:,5]
altitude = np.sqrt(xout**2 + yout**2 + zout**2) -Rplanet
velout = np.sqrt(velxout**2 + velyout**2 + velzout**2)


###Air density / Altitude
test_altitude = np.linspace(0, 600000, 100)
test_rho = aeroModel.getDensity(test_altitude)
plt.figure()
plt.plot(test_altitude, test_rho, 'b-')
plt.xlabel('Air Density (Kg/m^3)')
plt.ylabel('Altitude')
plt.grid()

###Altitude
plt.figure()
plt.plot(tout, altitude)
plt.xlabel('Time (s)')
plt.ylabel('Altitude (m)')
plt.grid()

###Velocity
plt.figure()
plt.plot(tout, velout)
plt.xlabel('Time (s)')
plt.ylabel('Total Speed (m/s)')
plt.grid()

####3D Orbit
u = np.linspace(0, 2*np.pi, 30) #longitude
v = np.linspace(0, np.pi, 30) #latitude
xsphere = Rplanet*np.outer(np.cos(u), np.sin(v))
ysphere = Rplanet*np.outer(np.sin(u), np.sin(v))
zsphere = Rplanet*np.outer(np.ones(np.size(u)), np.cos(v))

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.set_title('3D Cubesat Orbit Simulation')

ax.plot_surface(xsphere, ysphere, zsphere, color='b', alpha=0.3) #Planet

ax.plot(xout, yout, zout, 'r', label='Orbit')#Orbit



plt.legend()
ax.set_aspect('equal')
plt.show()

生成的图表： #

解决方案： #

要精准捕获碰撞时刻并终止仿真，我们可以利用scipy.integrate.odeint的事件检测功能，定义一个碰撞触发的终止条件。具体实现如下：

1. 定义碰撞事件函数

这个函数会返回卫星到行星中心的距离与行星半径的差值，当该值从正变负时（即卫星接触表面），触发终止事件：

def collision_event(state, t):
    x, y, z = state[0], state[1], state[2]
    r = np.sqrt(x**2 + y**2 + z**2)
    # 返回值为0时触发事件，direction=-1表示从正到负的跨越
    return r - Rplanet

# 设置事件属性：触发时立即终止积分
collision_event.terminal = True
collision_event.direction = -1

2. 修改仿真运行逻辑

调用odeint时传入事件参数，获取碰撞发生的时间和状态，并过滤掉碰撞后的无效数据：

# 运行积分并检测碰撞事件
stateout, info = sci.odeint(Derivatives, stateInitial, tout, full_output=True, events=collision_event)

# 处理碰撞结果
if info['t_events'][0].size > 0:
    collision_time = info['t_events'][0][0]
    print(f"碰撞发生在时间：{collision_time:.2f}秒")
    # 只保留碰撞前的有效数据
    mask = tout <= collision_time
    tout = tout[mask]
    stateout = stateout[mask]
else:
    print("在仿真时间窗口内未发生碰撞")

3. 优化辅助逻辑（可选）

在空气动力学计算中加入判断，避免碰撞后进行无效的密度计算：

###Aerodynamics
altitude = r - Rplanet
if altitude < 0:
    rho = 0
else:
    rho = aeroModel.getDensity(altitude)

4. 完整修改后的代码

import numpy as np
import scipy.integrate as sci
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

plt.close()

G = 6.6742e-11

###Planets
###Earth
name = 'Earth'
Rplanet= 6357000 #m
mplanet= 5.972e24 #Kg

###Cube satellite
mass = 1 #Kg
S = 1/100 #Transversal area
CD = 2

###Initial Conditions
x0 = 0
y0 = Rplanet + 350000
z0 = 0
r0 = Rplanet + 350000
velx0 = 0
vely0 = 0
velz0 = np.sqrt((G*mplanet)/r0)
period = 1e5

###AeroDynamics Model
class Aerodynamics():
    def __init__(self, name):
        self.name = name
        if name == 'Earth':
            self.beta = 0.1354/1000 #Density constant
            self.rhos = 1.225 #kg/m^3
    def getDensity(self, altitude):
        if self.name == 'Earth':
            rho = self.rhos*np.exp(-altitude*self.beta)
            return rho

aeroModel = Aerodynamics(name)

###Gravity model
def gravity(x, y, z):
    global Rplanet, mplanet, G

    r = np.sqrt(x**2 + y**2 + z**2)

    if r < Rplanet:
        accelx = 0
        accely = 0
        accelz = 0
    else:
        accelx = (G*mplanet)*x/(r**3)
        accely = (G*mplanet)*y/(r**3)
        accelz = (G*mplanet)*z/(r**3)
    return np.asarray([accelx, accely, accelz]), r

###Differential equation of the system
def Derivatives(state, t):
    global mass

    x = state[0]
    y = state[1]
    z = state[2]
    velx = state[3]
    vely = state[4]
    velz = state[5]

    xdot = velx
    ydot = vely
    zdot = velz

    ###Forces
    ###Gravity
    accel, r = gravity(x, y, z)
    GravityF = -accel*mass

    ###Aerodynamics
    altitude = r - Rplanet
    if altitude < 0:
        rho = 0
    else:
        rho = aeroModel.getDensity(altitude)
    V = np.sqrt(velx**2 + vely**2 + velz**2)
    qinf = -1/2*rho*abs(V)*S
    aeroF = qinf*CD*np.asarray([velx, vely, velz])
    
    Forces = GravityF + aeroF

    ddot = Forces/mass

    statedot = np.asarray([xdot, ydot, zdot, ddot[0], ddot[1], ddot[2]])
    return statedot

###碰撞事件检测函数
def collision_event(state, t):
    x, y, z = state[0], state[1], state[2]
    r = np.sqrt(x**2 + y**2 + z**2)
    return r - Rplanet

collision_event.terminal = True
collision_event.direction = -1

###Run simulation
stateInitial = np.asarray([x0, y0, z0, velx0, vely0, velz0])
tout = np.linspace(0, period, 1000)

stateout, info = sci.odeint(Derivatives, stateInitial, tout, full_output=True, events=collision_event)

# 处理碰撞结果
collision_time = None
if info['t_events'][0].size > 0:
    collision_time = info['t_events'][0][0]
    print(f"碰撞发生在时间：{collision_time:.2f}秒")
    mask = tout <= collision_time
    tout = tout[mask]
    stateout = stateout[mask]
else:
    print("在仿真时间窗口内未发生碰撞")

# 提取结果数据
xout = stateout[:,0]
yout = stateout[:,1]
zout = stateout[:,2]
velxout = stateout[:,3]
velyout = stateout[:,4]
velzout = stateout[:,5]
altitude = np.sqrt(xout**2 + yout**2 + zout**2) - Rplanet
velout = np.sqrt(velxout**2 + velyout**2 + velzout**2)

###Air density / Altitude plot
test_altitude = np.linspace(0, 600000, 100)
test_rho = aeroModel.getDensity(test_altitude)
plt.figure()
plt.plot(test_altitude, test_rho, 'b-')
plt.xlabel('Altitude (m)')
plt.ylabel('Air Density (Kg/m^3)')
plt.grid()

###Altitude plot
plt.figure()
plt.plot(tout, altitude)
plt.xlabel('Time (s)')
plt.ylabel('Altitude (m)')
plt.grid()
if collision_time:
    plt.axvline(x=collision_time, color='r', linestyle='--', label=f'碰撞时刻: {collision_time:.2f}s')
    plt.legend()

###Velocity plot
plt.figure()
plt.plot(tout, velout)
plt.xlabel('Time (s)')
plt.ylabel('Total Speed (m/s)')
plt.grid()
if collision_time:
    plt.axvline(x=collision_time, color='r', linestyle='--', label=f'碰撞时刻: {collision_time:.2f}s')
    plt.legend()

####3D Orbit plot
u = np.linspace(0, 2*np.pi, 30) #longitude
v = np.linspace(0, np.pi, 30) #latitude
xsphere = Rplanet*np.outer(np.cos(u), np.sin(v))
ysphere = Rplanet*np.outer(np.sin(u), np.sin(v))
zsphere = Rplanet*np.outer(np.ones(np.size(u)), np.cos(v))

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.set_title('3D Cubesat Orbit Simulation')

ax.plot_surface(xsphere, ysphere, zsphere, color='b', alpha=0.3) #Planet
ax.plot(xout, yout, zout, 'r', label='Orbit')#Orbit

plt.legend()
ax.set_aspect('equal')
plt.show()

修改后的代码会在卫星碰撞到行星表面时自动终止仿真，你可以直接获取精准的碰撞时刻，所有图表也只会显示到碰撞发生前的有效数据，不会再出现卫星进入行星内部的无效仿真结果。

备注：内容来源于stack exchange，提问作者Diniz Vitor