Завершение моделирования столкновения с поверхностью ⇐ Python

[Anonymous]

Я создаю симуляцию затухания орбиты и не знаю, как остановить/прервать симуляцию при столкновении объекта с поверхностью. Он продолжает моделировать поверхность планеты и не дает мне точного момента столкновения с поверхностью.
Вот код:

Код: Выделить всё

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

###Panets
###Earth
name = 'Earth'
Rplanet= 6357000 #K
mplanet= 5.972e24 #Kg

###Cube satelite
mass = 1 #Kg
S = 1/100 #Trasnversal 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 #Desity 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) #logitude
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()

И графики:







Подробнее здесь: https://stackoverflow.com/questions/799 ... -collision