ARTICLE AD BOX
I'm Making an orbit decay simulation and i dont know how can i make the simulation stop/break when the object collides with the surface. It keeps simulating inside the planet's surface and doesnt give me the precise moment of colision to the surface.
Here is the Code:
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()And the Plots:
