from skyfield.api import load import numpy as np from matplotlib import pyplot as plt from skyfield.constants import AU_KM from skyfield.vectorlib import VectorFunction from spktype21 import SPKType21 #Mapping SunID = 10 MercuryID = 0 VenusID = 1 EarthID = 2 MarsID = 3 JupiterID = 4 SaturnID = 5 UranusID = 6 NeptuneID = 7 PlutoID = 8 OumuamuaID = 9 MDID = 11 #Parameters dt = 1 tsteps = 80 N = 12 day0 = 14 #Start of observations day1 = 94 #One day past observations year0 = 2017 GM = 4.0*np.pi**2/(365.25)**2 froot = './Data/24_01_2025b_t' mN2 = 28 #Mass of N2 gas in grams/mole R = 8.314 #Joules/mole*Kelvin #km to AU conversion alpha = 149597870.691 #Masses #Unit is 10^24 kg m = np.ones(N, dtype=float) m[SunID] = 1988400 m[MercuryID] = 0.330 m[VenusID] = 4.87 m[EarthID] = 5.97 m[MarsID] = 0.642 m[JupiterID] = 1898 m[SaturnID] = 568 m[UranusID] = 86.8 m[NeptuneID] = 102 m[PlutoID] = 0.0130 m[MDID] = 1E-16 PlanetaryObjects = (MercuryID, VenusID, EarthID, MarsID, JupiterID, SaturnID, UranusID, NeptuneID, PlutoID) #Times # Load the ephemeris data ts = load.timescale() planets = load('./Ephemeris/de421.bsp') kernel = SPKType21.open('./Ephemeris/3788040_observed.bsp') times =[] for d in range(day0,day1): times.append(ts.utc(year0, 10, d)) #Other data structures #x, y, and z are positions of simulated `Oumuamua #vx, vy, and vz are velocities of simulated `Oumuamua #Ovx, Ovy, and Ovz are observed velocities of `Oumuamua #ax, ay, and az are the simulated accelerations of `Oumuamua due to gravity #aFix_x, aFix_y, aFix_z are corrected accelerations of `Oumuamua based on observational data x = np.zeros(tsteps*N,dtype=np.float64).reshape(N,tsteps) y = np.zeros(tsteps*N,dtype=np.float64).reshape(N,tsteps) z = np.zeros(tsteps*N,dtype=np.float64).reshape(N,tsteps) vx = np.zeros(tsteps,dtype=np.float64) vy = np.zeros(tsteps,dtype=np.float64) vz = np.zeros(tsteps,dtype=np.float64) Ovx = np.zeros(tsteps,dtype=np.float64) Ovy = np.zeros(tsteps,dtype=np.float64) Ovz = np.zeros(tsteps,dtype=np.float64) aFix_x = np.zeros(tsteps,dtype=np.float64) aFix_y = np.zeros(tsteps,dtype=np.float64) aFix_z = np.zeros(tsteps,dtype=np.float64) diffV = np.zeros(tsteps, dtype=np.float64) MassRatio = np.zeros(tsteps, dtype=np.float64) MassRatio20 = np.zeros(tsteps, dtype=np.float64) MassRatio40 = np.zeros(tsteps, dtype=np.float64) MassRatio30 = np.zeros(tsteps, dtype=np.float64) ax = 0 ay = 0 az = 0 #Computing `Oumuamua's initial position and velocity p, v = kernel.compute_type21(10, 3788040, times[0].tt) x[MDID, 0] = p[0]/alpha y[MDID, 0] = p[1]/alpha z[MDID, 0] = p[2]/alpha vx[0] = v[0]*24*60*60/alpha vy[0] = v[1]*24*60*60/alpha vz[0] = v[2]*24*60*60/alpha def Dist(x, y, z): return np.sqrt(x**2+y**2+z**2) def SunGrav(x,y,z): r = Dist(x,y,z) a = -(GM)/r**3 return x*a, y*a, z*a def AstroGrav(x,y,z,ax,ay,az,m,t): ax, ay, az = SunGrav(x[MDID,t],y[MDID,t],z[MDID,t]) for j in PlanetaryObjects: dx = x[MDID,t]-x[j,t] dy = y[MDID,t]-y[j,t] dz = z[MDID,t]-z[j,t] r = Dist(dx, dy, dz) if r < 0.0001: print(j) a = -(GM*(m[j]/m[SunID]))/r**3 ax = a*dx+ax ay = a*dy+ay az = a*dz+az return ax, ay, az def PositionCalculator(times,x,y,z): sun = planets['sun'] mercury = planets['mercury barycenter'] venus = planets['venus barycenter'] earth = planets['earth barycenter'] mars = planets['mars barycenter'] jupiter = planets['jupiter barycenter'] saturn = planets['saturn barycenter'] uranus = planets['uranus barycenter'] neptune = planets['neptune barycenter'] pluto = planets['pluto barycenter'] mercury_orbits = [] venus_orbits = [] earth_orbits = [] mars_orbits = [] jupiter_orbits = [] saturn_orbits = [] uranus_orbits = [] neptune_orbits = [] pluto_orbits = [] for t in times: sun_pos = sun.at(t).position.au mercury_pos = mercury.at(t).position.au venus_pos = venus.at(t).position.au earth_pos = earth.at(t).position.au mars_pos = mars.at(t).position.au jupiter_pos = jupiter.at(t).position.au saturn_pos = saturn.at(t).position.au uranus_pos = uranus.at(t).position.au neptune_pos = neptune.at(t).position.au pluto_pos = pluto.at(t).position.au mercury_orbits.append(mercury_pos - sun_pos) venus_orbits.append(venus_pos - sun_pos) earth_orbits.append(earth_pos - sun_pos) mars_orbits.append(mars_pos - sun_pos) jupiter_orbits.append(jupiter_pos - sun_pos) saturn_orbits.append(saturn_pos - sun_pos) uranus_orbits.append(uranus_pos - sun_pos) neptune_orbits.append(neptune_pos - sun_pos) pluto_orbits.append(pluto_pos - sun_pos) Nt = len(times) for tslice in range(0, Nt): x[MercuryID, tslice] = mercury_orbits[tslice][0] y[MercuryID, tslice] = mercury_orbits[tslice][1] z[MercuryID, tslice] = mercury_orbits[tslice][2] x[VenusID, tslice] = venus_orbits[tslice][0] y[VenusID, tslice] = venus_orbits[tslice][1] z[VenusID, tslice] = venus_orbits[tslice][2] x[EarthID, tslice] = earth_orbits[tslice][0] y[EarthID, tslice] = earth_orbits[tslice][1] z[EarthID, tslice] = earth_orbits[tslice][2] x[MarsID, tslice] = mars_orbits[tslice][0] y[MarsID, tslice] = mars_orbits[tslice][1] z[MarsID, tslice] = mars_orbits[tslice][2] x[JupiterID, tslice] = jupiter_orbits[tslice][0] y[JupiterID, tslice] = jupiter_orbits[tslice][1] z[JupiterID, tslice] = jupiter_orbits[tslice][2] x[SaturnID, tslice] = saturn_orbits[tslice][0] y[SaturnID, tslice] = saturn_orbits[tslice][1] z[SaturnID, tslice] = saturn_orbits[tslice][2] x[UranusID, tslice] = uranus_orbits[tslice][0] y[UranusID, tslice] = uranus_orbits[tslice][1] z[UranusID, tslice] = uranus_orbits[tslice][2] x[NeptuneID, tslice] = neptune_orbits[tslice][0] y[NeptuneID, tslice] = neptune_orbits[tslice][1] z[NeptuneID, tslice] = neptune_orbits[tslice][2] x[PlutoID, tslice] = pluto_orbits[tslice][0] y[PlutoID, tslice] = pluto_orbits[tslice][1] z[PlutoID, tslice] = pluto_orbits[tslice][2] def GetOumuamua(times): oumuamua_movement = [] for t in times: oumuamua_pos, junk = kernel.compute_type21(10, 3788040, t.tt) oumuamua_movement.append(oumuamua_pos/alpha) Nt = len(times) for tslice in range(0, Nt): x[OumuamuaID, tslice] = oumuamua_movement[tslice][0] y[OumuamuaID, tslice] = oumuamua_movement[tslice][1] z[OumuamuaID, tslice] = oumuamua_movement[tslice][2] def GetOumuamuaV(times): oumuamua_v = [] for t in times: oumuamua_pos, v = kernel.compute_type21(10, 3788040, t.tt) oumuamua_v.append(v*24*60*60/alpha) Nt = len(times) for tslice in range(0, Nt): Ovx[tslice] = oumuamua_v[tslice][0] Ovy[tslice] = oumuamua_v[tslice][1] Ovz[tslice] = oumuamua_v[tslice][2] def DumpCoords(times, froot, x, y, z): Nt = len(times) for tslice in range(0,Nt): fname = froot+str(tslice)+'.csv' with open(fname, "w") as f: f.write('planet,x, y, z\n') f.write('10,0,0,0\n') for object in PlanetaryObjects: f.write(str(object) + ','+ str(x[object,tslice]) + ',' + str(y[object,tslice]) + ',' + str(z[object,tslice])+'\n') f.write(str(OumuamuaID) + ','+ str(x[OumuamuaID,tslice]) + ',' + str(y[OumuamuaID,tslice]) + ',' + str(z[OumuamuaID,tslice])+'\n') f.write(str(MDID) + ','+ str(x[MDID,tslice]) + ',' + str(y[MDID,tslice]) + ',' + str(z[MDID,tslice])+'\n') # Calculates the positions for all the planets PositionCalculator(times,x,y,z) #Loads Oumuamua observed data and calculates its positions GetOumuamua(times) #Calculates the effect of the planets on Oumuamua at the first timestep ax, ay, az = AstroGrav(x,y,z,ax,ay,az,m,0) print(ax,ay,az) #fvx, fvy, and fvz are the corrected velocities of `Oumuamua based on observational data fvx = np.zeros(tsteps,dtype=np.float64) fvy = np.zeros(tsteps,dtype=np.float64) fvz = np.zeros(tsteps,dtype=np.float64) fvx[0] = v[0]*24*60*60/alpha fvy[0] = v[1]*24*60*60/alpha fvz[0] = v[2]*24*60*60/alpha def Sim():#Fix): #aFix_x, aFix_y, aFix_z = Fix ax = 0 ay = 0 az = 0 ax, ay, az = AstroGrav(x,y,z,ax,ay,az,m,0) aFix_x[0] = 1/dt*( (Ovx[1]-vx[1])) aFix_y[0] = 1/dt*( (Ovy[1]-vy[1])) aFix_z[0] = 1/dt*( (Ovz[1]-vz[1])) ax = ax+ aFix_x[0] ay = ay+ aFix_y[0] az = az+ aFix_z[0] for t in range(1,tsteps): x[MDID,t] = x[MDID,t-1] + fvx[t-1]*dt + 0.5*ax*dt*dt y[MDID,t] = y[MDID,t-1] + fvy[t-1]*dt + 0.5*ay*dt*dt z[MDID,t] = z[MDID,t-1] + fvz[t-1]*dt + 0.5*az*dt*dt fvx[t] = fvx[t-1] + 0.5*dt*ax fvy[t] = fvy[t-1] + 0.5*dt*ay fvz[t] = fvz[t-1] + 0.5*dt*az ax, ay, az = AstroGrav(x,y,z,ax,ay,az,m,t) fvx[t] = fvx[t] + 0.5*dt*ax fvy[t] = fvy[t] + 0.5*dt*ay fvz[t] = fvz[t] + 0.5*dt*az aFix_x[t] = 1/dt*( (Ovx[t]-fvx[t])) aFix_y[t] = 1/dt*( (Ovy[t]-fvy[t])) aFix_z[t] = 1/dt*( (Ovz[t]-fvz[t])) ax = ax+ aFix_x[t] ay = ay+ aFix_y[t] az = az+ aFix_z[t] fvx[t] = fvx[t] + 1*dt*aFix_x[t] fvy[t] = fvy[t] + 1*dt*aFix_y[t] fvz[t] = fvz[t] + 1*dt*aFix_z[t] #Simulation of `Oumuamua positions, velocities and accelerations due to gravity # Calculates the positions for all the planets PositionCalculator(times,x,y,z) #Loads Oumuamua observed data and calculates its positions GetOumuamua(times) #Calculates the effect of the planets on Oumuamua at the first timestep ax, ay, az = AstroGrav(x,y,z,ax,ay,az,m,0) for t in range(1,tsteps): x[MDID,t] = x[MDID,t-1] + vx[t-1]*dt + 0.5*ax*dt*dt y[MDID,t] = y[MDID,t-1] + vy[t-1]*dt + 0.5*ay*dt*dt z[MDID,t] = z[MDID,t-1] + vz[t-1]*dt + 0.5*az*dt*dt vx[t] = vx[t-1] + 0.5*dt*ax vy[t] = vy[t-1] + 0.5*dt*ay vz[t] = vz[t-1] + 0.5*dt*az ax, ay, az = AstroGrav(x,y,z,ax,ay,az,m,t) vx[t] = vx[t] + 0.5*dt*ax vy[t] = vy[t] + 0.5*dt*ay vz[t] = vz[t] + 0.5*dt*az DumpCoords(times, froot, x, y, z) #Runs simulation where simulated acceleration is corrected at each time step GetOumuamuaV(times) Sim() Err2 = 0.0 for t in range(0,tsteps): Vtmp = Dist(fvx[t],fvy[t],fvz[t]) OVtmp = Dist(Ovx[t],Ovy[t],Ovz[t]) Err2 = Err2+((Vtmp-OVtmp)*alpha)**2 Err2