pythonphysics

Why is my Nbody simulator not printing out orbital times past 3 bodies?


Here is my code for simulating planetary orbits. When I have the bodies list set up with only the Earth, Sun and Jupiter, my code works well and prints out a reasonably accurate time for Jupiter and Earth's orbit. However, when I add Saturn into the bodies list, I get a value of 43200 seconds for both Jupiter and Saturn's orbit. What's strange though is that Earth's orbit is printing out accurately even with Saturn in the list of bodies.

In addition to this, the plot I'm getting is correct, as if I run the simulation for, say, 28 years, Saturn shows as not having a completed orbit. If I run the simulation for 29.5 years, Saturn just about completes its orbit, which is accurate. So the gravitational calculation is still working. Perhaps the index is failing somehow when another body is added? But I'm not sure why.

It seems the line

simulation.run(30*u.year,6*u.hr)

is causing this. Whenever I change the time-step to around 24*u.hr, it prints out accurate adjacent times for only Jupiter, but not accurate enough. When I change the time-step to 10 days, it prints out a value of 10600 for Saturn's orbital period. This isn't accurate enough however.

import numpy as np 
import matplotlib.pyplot as plt 
import astropy.units as u 
import astropy.constants as c 
import sys 
import time
from mpl_toolkits.mplot3d import Axes3D

#making a class for Celestial Objects
class CelestialObjects():
    def __init__(self,mass,pos_vec,vel_vec,name=None, has_units=True):
        self.name=name
        self.has_units=has_units
        if self.has_units:
            self.mass=mass.cgs
            self.pos=pos_vec.cgs.value
            self.vel=vel_vec.cgs.value
        else:
            self.mass=mass 
            #3d array for position of body in 3d space in AU
            self.pos=pos_vec 
            #3d array for velocity of body in 3d space in km/s
            self.vel=vel_vec
        
    def return_vec(self):
        return np.concatenate((self.pos,self.vel))
    def return_name(self):
        return self.name
    def return_mass(self):
        if self.has_units:
            return self.mass.cgs.value
        else:
            return self.mass

v_earth=(((c.G*1.98892E30)/1.495978707E11)**0.5)/1000
v_jupiter=(((c.G*1.98892E30)/7.779089276E11)**0.5)/1000
v_saturn=(((c.G*1.98892E30)/1.421179772E12)**0.5)/1000

#set up first instance of a celestial object, Earth
Earth=CelestialObjects(name='Earth',
                       pos_vec=np.array([0,1,0])*u.AU,
                       vel_vec=np.array([v_earth.value,0,0])*u.km/u.s,
                       mass=1.0*c.M_earth)
#set up second instance of a celestial object, the Sun
Sun=CelestialObjects(name='Sun',
                     pos_vec=np.array([0,0,0])*u.AU,
                     vel_vec=np.array([0,0,0])*u.km/u.s,
                     mass=1*u.Msun)
Jupiter=CelestialObjects(name='Jupiter', 
                         pos_vec=np.array([0,5.2,0])*u.AU, 
                         vel_vec=np.array([v_jupiter.value,0,0])*u.km/u.s,
                         mass=317.8*c.M_earth)
Saturn=CelestialObjects(name='Saturn',
                        pos_vec=np.array([0,9.5,0])*u.AU, 
                        vel_vec=np.array([v_saturn.value,0,0])*u.km/u.s,
                        mass=95.16*c.M_earth)
                       
bodies=[Earth,Sun,Jupiter,Saturn]
#making a class for system
class Simulation():
    def __init__(self,bodies,has_units=True):
        self.has_units=has_units
        self.bodies=bodies
        self.Nbodies=len(self.bodies)
        self.Ndim=6
        self.quant_vec=np.concatenate(np.array([i.return_vec() for i in self.bodies]))
        self.mass_vec=np.array([i.return_mass() for i in self.bodies])
        self.name_vec=[i.return_name() for i in self.bodies]
        
    def set_diff_eqs(self,calc_diff_eqs,**kwargs):
        self.diff_eqs_kwargs=kwargs
        self.calc_diff_eqs=calc_diff_eqs
        
    
    def rk4(self,t,dt):
        k1= dt* self.calc_diff_eqs(t,self.quant_vec,self.mass_vec,**self.diff_eqs_kwargs)
        k2=dt*self.calc_diff_eqs(t+dt*0.5,self.quant_vec+0.5*k1,self.mass_vec,**self.diff_eqs_kwargs)
        k3=dt*self.calc_diff_eqs(t+dt*0.5,self.quant_vec+0.5*k2,self.mass_vec,**self.diff_eqs_kwargs)
        k4=dt*self.calc_diff_eqs(t+dt,self.quant_vec+k3,self.mass_vec,**self.diff_eqs_kwargs)
        
        y_new=self.quant_vec+((k1+2*k2+2*k3+k4)/6)
        return y_new
    
    def run(self,T,dt,t0=0):
        if not hasattr(self,'calc_diff_eqs'):
            raise AttributeError('You must set a diff eq solver first.')
        if self.has_units:
            try:
                _=t0.unit
            except:
                t0=(t0*T.unit).cgs.value
            T=T.cgs.value
            dt=dt.cgs.value
        
        self.history=[self.quant_vec]
        clock_time=t0
        nsteps=int((T-t0)/dt)
        start_time=time.time()
        orbit_completed=False
        orbit_start_time=clock_time
        initial_position=self.quant_vec[:3]
        min_distance=float('inf')
        min_distance_time=0
        count=0
        
        for step in range(nsteps):
            sys.stdout.flush()
            sys.stdout.write('Integrating: step = {}/{}| Simulation Time = {}'.format(step,nsteps,round(clock_time,3))+'\r')
            y_new=self.rk4(0,dt)
            self.history.append(y_new)
            self.quant_vec=y_new
            clock_time+=dt
            #checking if planet has completed an orbit
            current_position=self.quant_vec[:3] #**THIS IS WHERE ORBIT TIME CALCULATED** To explain, this is earths position, the first three elements of the vector, the next three elements are its velocity, and then the next three are the suns position vectors and so on. Saturns index would be self.quant_vec[18:21]. 

            distance_to_initial=np.linalg.norm(current_position-initial_position)
            if distance_to_initial<min_distance and orbit_completed is False:
                min_distance=distance_to_initial
                min_distance_time=clock_time
                if count==1:
                    orbit_completed=True
                count+=1
                
                
                
                
            
        runtime=time.time()-start_time
        print(clock_time)
        print('\n')
        print('Simulation completed in {} seconds'.format(runtime))
        print(f'Minimum distance reached at time: {min_distance_time:.3f} seconds. Minimum distance: {min_distance:.3e}')
        self.history=np.array(self.history)
    
def nbody_solver(t,y,masses):
    N_bodies=int(len(y)/6)
    solved_vector=np.zeros(y.size)
    distance=[]
    for i in range(N_bodies):
        ioffset=i * 6
        for j in range(N_bodies):
            joffset=j * 6
            solved_vector[ioffset]=y[ioffset+3]
            solved_vector[ioffset+1]=y[ioffset+4]
            solved_vector[ioffset+2]=y[ioffset+5]
            if i != j:
                dx= y[ioffset]-y[joffset]
                dy=y[ioffset+1]-y[joffset+1]
                dz=y[ioffset+2]-y[joffset+2]
                r=(dx**2+dy**2+dz**2)**0.5
                ax=(-c.G.cgs*masses[j]/r**3)*dx
                ay=(-c.G.cgs*masses[j]/r**3)*dy
                az=(-c.G.cgs*masses[j]/r**3)*dz
                ax=ax.value
                ay=ay.value
                az=az.value
                solved_vector[ioffset+3]+=ax
                solved_vector[ioffset+4]+=ay
                solved_vector[ioffset+5]+=az
    return solved_vector

simulation=Simulation(bodies)
simulation.set_diff_eqs(nbody_solver)
simulation.run(30*u.year,12*u.hr)


earth_position = simulation.history[:, :3]  # Extracting position for Earth
sun_position = simulation.history[:, 6:9]    # Extracting position for Sun
jupiter_position=simulation.history[:, 12:15] #Jupiter position
saturn_position=simulation.history[:, 18:21] #Saturn position

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the trajectories
ax.plot(earth_position[:, 0], earth_position[:, 1], earth_position[:, 2], label='Earth')
ax.plot(sun_position[:, 0], sun_position[:, 1], sun_position[:, 2], label='Sun')
ax.plot(jupiter_position[:, 0], jupiter_position[:, 1], jupiter_position[:, 2], label='Jupiter')
ax.plot(saturn_position[:, 0], saturn_position[:, 1], saturn_position[:, 2], label='Saturn')

# Add labels and title
ax.set_xlabel('X (AU)')
ax.set_ylabel('Y (AU)')
ax.set_zlabel('Z (AU)')
ax.set_title('Trajectories of Earth and Jupiter Around the Sun')
ax.scatter([0], [0], [0], marker='o', color='yellow', s=200, label='Sun')  # Marking the Sun at the origin
ax.scatter(earth_position[0, 0], earth_position[0, 1], earth_position[0, 2], marker='o', color='blue', s=50, label='Earth')  # Marking the initial position of Earth
ax.scatter(jupiter_position[0, 0], jupiter_position[0, 1], jupiter_position[0, 2], marker='o', color='green', s=100, label='Jupiter')  # Marking the initial position of Earth
ax.scatter(saturn_position[0, 0], saturn_position[0, 1], saturn_position[0, 2], marker='o', color='red', s=80, label='Saturn')  # Marking the initial position of Earth
# Add a legend
ax.legend()

# Show the plot
plt.show()

Solution

  • The reason for this behaviour is that the if distance_to_initial<min_distance condition for the time logging is not triggered correctly in this case.

    The specific cause of this is the "min_distance", which might get "hopped over" in certain circumstances (e.g. due to increasing step distance), and then only the first min_distance_time value is used, without ever "completing the orbit".

    To fix this, instead of using a min_distance, i would instead check if the initial position is in between the previous and the current position.

    Here's a run method which implements this approach:

        def run(self, T, dt, t0=0):
            if not hasattr(self, "calc_diff_eqs"):
                raise AttributeError("You must set a diff eq solver first.")
            if self.has_units:
                try:
                    _ = t0.unit
                except:
                    t0 = (t0 * T.unit).cgs.value
                T = T.cgs.value
                dt = dt.cgs.value
    
            self.history = [self.quant_vec]
            clock_time = t0
            nsteps = int((T - t0) / dt)
            start_time = time.time()
    
            initial_position = self.quant_vec[18:21]
    
            min_distance = None
            min_distance_time = None
    
            for step in range(nsteps):
                sys.stdout.flush()
                sys.stdout.write(
                    "Integrating: step = {}/{}| Simulation Time = {}".format(
                        step, nsteps, round(clock_time, 3)
                    )
                    + "\r"
                )
                y_new = self.rk4(0, dt)
                y_old = self.quant_vec
                self.history.append(y_new)
                self.quant_vec = y_new
                clock_time += dt
                current_position = self.quant_vec[18:21]
                last_position = y_old[18:21]
    
                if step == 0 or min_distance_time is not None:
                    # do not check orbit criteria for initial step
                    # or if orbit already completed
                    continue  
    
                # checking if planet has completed an orbit, i.e. if origin is between current and last position
                last_to_initial = np.linalg.norm(last_position - initial_position)
                current_to_initial = np.linalg.norm(current_position - initial_position)
                current_to_last = np.linalg.norm(current_position - last_position)
                has_passed_initial_position = (
                    last_to_initial <= current_to_last
                    and current_to_initial <= current_to_last
                )
                if has_passed_initial_position:
                    min_distance = current_to_initial
                    min_distance_time = clock_time
    
            runtime = time.time() - start_time
            print("\n")
            print("Simulation completed in {} seconds".format(runtime))
            print(
                f"Minimum distance reached at time: {min_distance_time:.3f} seconds. Minimum distance: {min_distance:.3e}"
            )
            self.history = np.array(self.history)