<p>I am relatively new to astronomical programming and I am trying to code a <strong>simulation of the solar system</strong> as simplified as possible on python. To do so, I am using numpy in order to have one array containing in each row the x,y,z coordinates of position of each celestial body. So if N=5 bodies, the positions matrix will be sized (5,3). The initial positions and velocities of the bodies are uploading with the <strong>Spiceypy</strong> library (SPICE ephemeris). Then, I'm integrating using the <strong>leapfrog</strong> method. However, my code seems to run well until I plot the resulting positions: I keep on obtaining <strong>straight lines</strong> instead of elliptical orbits...</p>
<p><a href="https://i.sstatic.net/9s2e6.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9s2e6.png" alt="enter image description here" /></a></p>
<p><strong>I am definitely doing something wrong, but I don't know what. Do you know what I am missing here? I suspect the problem is in the leapfrog integration.</strong></p>
<p>Thank you for any kind of help you can provide!</p>
<p>Here is my code:</p>
<pre><code>    # -*- coding: utf-8 -*-
&quot;&quot;&quot;
Created on Tue Apr 26 12:35:05 2022

@author: Marwa Kadhem
&quot;&quot;&quot;


import numpy as np
import matplotlib.pyplot as plt
import spiceypy as spice


spice.tkvrsn('TOOLKIT')
spice.furnsh('de440.bsp')
spice.furnsh('p10-a.bsp')
spice.furnsh('naif0012.tls')



def Ephemeris(body):
    
    #-------------------------------------------------------------------------
    # INPUT
    # body : name of the body (type string)
    #-------------------------------------------------------------------------
    # OUPUT
    # (pos,vel) : a tuple where index 0 is x,y,z of position (np array) and 
    #             index 1 is x,y,z of velocity (np array) of the body obtained 
    #             at time ref_time from NASA SPICE files
    #-------------------------------------------------------------------------
    
    ref_time = spice.spiceypy.str2et('2022 May 01 00:00:00.0000')
    
    pos_list = spice.spiceypy.spkezr(body,ref_time,'J2000','NONE','SUN')[0][0:3]
    pos_matrix = np.array([ [pos_list[0], pos_list[1],pos_list[2]] ])
    
    vel_list = spice.spiceypy.spkezr(body,ref_time,'J2000','NONE','SUN')[0][3:6]
    vel_matrix = np.array([ [vel_list[0], vel_list[1],vel_list[2]] ])
    
    return (pos_matrix,vel_matrix)


def NbodyProblem(pos, mass, G):
    
    #-----------------------------------------------------------------
    # INPUT
    # pos  : positions matrix size (N,3) 
    #        where each row i corresponds to the x,y,z coordinate of 
    #        one celestial body of the Solar System.
    # mass : masses matrix  (1,N)
    #        each row is the mass (kg) of one celestial body of the
    #        Solar System.
    # G    : Newton's Gravitational constant
    #-----------------------------------------------------------------
    # OUTPUT
    # acc  : accelerations matrix (N,3)
    #        where each row i corresponds to the ax,ay,az  
    #        accelerations of one celestial body of the Solar System.
    #-----------------------------------------------------------------
    
    x = [row[0] for row in pos] # x element of position matrix
    y = [row[1] for row in pos] # y element of position matrix
    z = [row[2] for row in pos] # x element of position matrix
    acc = np.empty((0,3)) 
    
    for j in range (0,len(x)):
        ax=0
        ay=0
        az=0
        
        for i in range (0,len(x)):
            if (i!=j):
                dx = x[i]-x[j]
                dy = y[i]-y[j]
                dz = z[i]-z[j]
                norm = np.sqrt(dx**2+dy**2+dz**2)
                ax -= (G*mass[i]*dx)/norm**3
                ay -= (G*mass[i]*dy)/norm**3
                az -= (G*mass[i]*dz)/norm**3
                
        a = np.array([ax,ay,az])
        acc = np.append(acc,[a],axis=0)
        
    return acc


def leapfrog(pos0, vel0, mass, dt):
    
    #-----------------------------------------------------------------
    # INPUT
    # pos0 : initial positions matrix size (N,3). 
    # vel0 : initial velocities matrix size (N,3).
    # mass : masses matrix (1,N).
    # time : time of the simulation ranging from start time (t_start)
    #        to end time (t_end).
    # dt   : timestep.
    #-----------------------------------------------------------------
    # OUTPUT
    # pos_register : dictionnary of all positions matrix size (N,3)
    #                found by leapfrog or Stormer-Verlet integration 
    #                method at each time of the simulation.
    #-----------------------------------------------------------------
    
    global G
    global t_start
    global t_end
    
    pos_register = {}
    
    pos = pos0
    vel = vel0
    acc = NbodyProblem(pos,mass,G)
    
    for t in np.arange(t_start, t_end, dt): # kick-drift-kick version

        vel = vel + acc * dt * 0.5       # kick
        pos = pos + vel * dt             # drift
        
        acc = NbodyProblem(pos,mass,G)   # update acceleration
        vel = vel + acc * dt * 0.5       # kick
        
        pos_register[t] = pos            # save in the register
        
    return pos_register




        
def Coordinates(register, body_id):
    
    #----------------------------------------------------------------
    # INPUT
    # register : Dictionary containing all the positions matrix 
    #            at all time t of the simulation.
    # body_id : identification number of the body, which corresponds 
    #           to the row number in the positions matrix. 
    #----------------------------------------------------------------
    # OUTPUT
    # (Lx,Ly,Lz) : 3 lists containig respectively all x, y and z
    #              values of the desired body (id_body).
    #----------------------------------------------------------------
    
    global t_start
    global t_end
    global dt
    
    Lx = []
    Ly = []
    Lz = []
    for t in np.arange(t_start, t_end, dt):
        matrix = register[t]
        Lx.append(matrix[body_id,0])
        Ly.append(matrix[body_id,1])
        Lz.append(matrix[body_id,2])
        
    return (Lx,Ly,Lz)    
#=======================================================================================================
#                                           Simulation parameters
#=======================================================================================================

N = 5                                                               # number of celestial bodies

t_start = spice.spiceypy.str2et('2022 May 01 00:00:00.0000')        # start time of the simulation

t_end = spice.spiceypy.str2et('2050 May 01 00:00:00.0000')          # end time of simulation

dt = 24*3600                                                           # timestep

G = 6.67430e-11*1e-9                                                    # Newton's Gravitational Constant

num_ts = int(np.ceil(t_end/dt))                                     # number of timesteps

#=======================================================================================================


# EPHEMERIS 2022 May 01 00:00:00.0000

# Sun
sun_id = 0
pos_sun = Ephemeris('SUN')[0]
vel_sun = Ephemeris('SUN')[1]
m_sun = 1.989e30


# Mercury 
mercury_id = 1
pos_mercury = Ephemeris('MERCURY')[0]
vel_mercury = Ephemeris('MERCURY')[1]
m_mercury = 0.330e24


# Mars
mars_id = 2
pos_mars = Ephemeris('MARS BARYCENTER')[0]  
vel_mars = Ephemeris('MARS BARYCENTER')[1]   
m_mars = 0.642e24


# Earth
earth_id = 3
pos_earth = Ephemeris('EARTH')[0]  
vel_earth = Ephemeris('EARTH')[1]  
m_earth = 5.97e24

# Venus
venus_id = 4
pos_venus = Ephemeris('VENUS')[0]    
vel_venus = Ephemeris('VENUS')[1]  
m_venus = 4.87e24


# Initial conditions mass, positions and velocities with ephemeris

M = np.array([m_sun,m_mercury,m_mars,m_earth,m_venus])

P = np.concatenate((pos_sun,pos_mercury,pos_mars,pos_earth,pos_venus),axis=0)

V = np.concatenate((vel_sun,vel_mercury,vel_mars,vel_earth,vel_venus),axis=0)

V -= np.mean(np.dot(M,V),0)/np.mean(M)


# Integration using the leapfrog integration (Störmer-Verlet method)

reg = leapfrog(P,V,M,dt)
print(reg)

# Generating lists containing x,y,z position values of each bodies
Lx_sun = Coordinates(reg,0)[0]
Ly_sun = Coordinates(reg,0)[1]
Lz_sun = Coordinates(reg,0)[2]
print(Lx_sun)

Lx_mercury = Coordinates(reg,1)[0]
Ly_mercury = Coordinates(reg,1)[1]
Lz_mercury = Coordinates(reg,1)[2]

Lx_mars = Coordinates(reg,2)[0]
Ly_mars = Coordinates(reg,2)[1]
Lz_mars = Coordinates(reg,2)[2]

Lx_earth = Coordinates(reg,3)[0]
Ly_earth = Coordinates(reg,3)[1]
Lz_earth = Coordinates(reg,3)[2]

Lx_venus = Coordinates(reg,4)[0]
Ly_venus = Coordinates(reg,4)[1]
Lz_venus = Coordinates(reg,4)[2]


# Ploting the results using Matplotlib

fig = plt.figure(&quot;Solar System simulation - Marwa Kadhem&quot;)

ax = fig.add_subplot(111, projection='3d')
ax.set_axis_off()
fig.set_facecolor('black')
ax.set_facecolor('black') 

ax.scatter(0,0,0, s=100, c=&quot;yellow&quot;,label='Sun')
ax.plot(Lx_mercury,Ly_mercury,Lz_mercury, label='Mercury')
ax.plot(Lx_mars,Ly_mars,Lz_mars, label='Mars')
ax.plot(Lx_earth,Ly_earth,Lz_earth, label='Earth')
ax.plot(Lx_venus,Ly_venus,Lz_venus, label='Venus')
ax.legend()

plt.show()
</code></pre>

CC BY-SA 4.0
2022-05-30T13:12:21.203Z
Hi everyone, I can't thank you enough for all your suggestions and help. I manage to solve the problem by following the steps you described and by taking into account all your replies.
mkvb_78
CC BY-SA 4.0
2022-05-30T13:15:32.643Z
If you want to check the result, it is fully available at https://github.com/MarwaKadhem/Solar-System-3D-simulation
mkvb_78
CC BY-SA 4.0
2022-05-30T13:16:49.567Z
And for the demo video: https://m.youtube.com/watch?v=g7ucAwWma2Y
mkvb_78
CC BY-SA 4.0
2022-05-02T11:57:10.533Z
First guess is to make sure all of your units are correct.  This sort of result is common when the force of gravity far exceeds the initial velocities.
Greg Miller
CC BY-SA 4.0
2022-05-02T16:01:37.927Z
In addition, numerical stability is much easier if you use AU-solar mass-years units, instead of SI units.  For example, in these units G=4pi.
James K
CC BY-SA 4.0
2022-05-02T16:45:40.100Z
finally, for testing purposes, it would be convenient to break the dependency on spice and  de440 (as these are large file dependencies.)  You are only using them to initialise the system.  So I'd suggest working out the initial values for may 2022 and hard coding them into the program.  When everything is working with the hard-coded values, you can easily switch back the de400 values.
James K
CC BY-SA 4.0
2022-05-02T17:06:13.207Z
Why do you have `G = 6.67430e-11*1e-9`. Where does that extra `*1e-9` come from? Everything else seems to be in mks units but this value. Unrelatedly, I'd suggest you post your code on something like [code review](https://codereview.stackexchange.com/). I think you could really benefit from someone reviewing your code from a coding perspective (aside from the science perspective you're asking for here).
zephyr
CC BY-SA 4.0
2022-05-03T00:05:00.173Z
Straight lines intersecting the Sun mean that gravity is too strong or speed is too low. A quick way to check is to calculate the *reduced* kinetic and potential energies $v^2/2$ and $-GM_{Sun}/r^2$ and see if the kinetic is about half of the absolute value of the potential energy since a bound orbit will have negative total energy. Next, if the lines radiate *away* from the Sun then you might be missing a minus sign somewhere. Also be extremely careful to never mix kilograms and seconds with kilometers, use only meters (perhaps the 1E-09 zephyr mentions comes from that?)
uhoh
CC BY-SA 4.0
2022-05-03T00:16:21.573Z
Finally, back up and first see if your method can solve a very familiar one-body problem. In a much shorter, simpler script, just put the Earth at $1.5 \times10^{11} \ \mathbf{\hat{x}}$ meters with a velocity of $2.9780 \times10^4 \ \mathbf{\hat{y}}$ meters per second, and confirm it goes around the $\mathbf{\hat{z}}$ axis  in roughly a circle and takes about a year to do so. All that said, here's a complicated n-body script that matches a JPL ephemeris pretty well https://space.stackexchange.com/a/23409/12102
uhoh
CC BY-SA 4.0
2022-05-03T00:27:19.787Z
fyi "one-body" means just calculate for the Earth and leave the Sun's central force at the origin. Once you get the integrator working and generating real, physical results, then work your way up slowly; add several bodies without interactions (n x 1-body) then finally add back all the $n(n-1$ two-body interactions i.e. n-body.
uhoh
CC BY-SA 4.0
2022-05-05T14:05:44.410Z
Since you're doing the calculations yourself, you might want to look at https://astronomy.stackexchange.com/questions/2416/where-can-i-find-a-set-of-data-of-the-initial-conditions-of-our-solar-system/13491#13491 and perhaps https://astronomy.stackexchange.com/questions/13488 (under "If you want to compute star/planet positions yourself, you have several options")
Barry Carter
<solar-system><orbital-mechanics><python><n-body-simulations>
Nbody problem solving on python with Leapfrog
mkvb_78
# Nbody problem solving on python with Leapfrog

(Asked by: username_0 on 2022-05-02)

I am relatively new to astronomical programming and I am trying to code a simulation of the solar system as simplified as possible on python. To do so, I am using numpy in order to have one array containing in each row the x,y,z coordinates of position of each celestial body. So if N=5 bodies, the positions matrix will be sized (5,3). The initial positions and velocities of the bodies are uploading with the Spiceypy library (SPICE ephemeris). Then, I'm integrating using the leapfrog method. However, my code seems to run well until I plot the resulting positions: I keep on obtaining straight lines instead of elliptical orbits...

I am definitely doing something wrong, but I don't know what. Do you know what I am missing here? I suspect the problem is in the leapfrog integration.
Thank you for any kind of help you can provide!
Here is my code:
    # -*- coding: utf-8 -*-
"""
Created on Tue Apr 26 12:35:05 2022

@author: Marwa Kadhem
"""


import numpy as np
import matplotlib.pyplot as plt
import spiceypy as spice


spice.tkvrsn('TOOLKIT')
spice.furnsh('de440.bsp')
spice.furnsh('p10-a.bsp')
spice.furnsh('naif0012.tls')



def Ephemeris(body):
    
    #-------------------------------------------------------------------------
    # INPUT
    # body : name of the body (type string)
    #-------------------------------------------------------------------------
    # OUPUT
    # (pos,vel) : a tuple where index 0 is x,y,z of position (np array) and 
    #             index 1 is x,y,z of velocity (np array) of the body obtained 
    #             at time ref_time from NASA SPICE files
    #-------------------------------------------------------------------------
    
    ref_time = spice.spiceypy.str2et('2022 May 01 00:00:00.0000')
    
    pos_list = spice.spiceypy.spkezr(body,ref_time,'J2000','NONE','SUN')[0][0:3]
    pos_matrix = np.array([ [pos_list[0], pos_list[1],pos_list[2]] ])
    
    vel_list = spice.spiceypy.spkezr(body,ref_time,'J2000','NONE','SUN')[0][3:6]
    vel_matrix = np.array([ [vel_list[0], vel_list[1],vel_list[2]] ])
    
    return (pos_matrix,vel_matrix)


def NbodyProblem(pos, mass, G):
    
    #-----------------------------------------------------------------
    # INPUT
    # pos  : positions matrix size (N,3) 
    #        where each row i corresponds to the x,y,z coordinate of 
    #        one celestial body of the Solar System.
    # mass : masses matrix  (1,N)
    #        each row is the mass (kg) of one celestial body of the
    #        Solar System.
    # G    : Newton's Gravitational constant
    #-----------------------------------------------------------------
    # OUTPUT
    # acc  : accelerations matrix (N,3)
    #        where each row i corresponds to the ax,ay,az  
    #        accelerations of one celestial body of the Solar System.
    #-----------------------------------------------------------------
    
    x = [row[0] for row in pos] # x element of position matrix
    y = [row[1] for row in pos] # y element of position matrix
    z = [row[2] for row in pos] # x element of position matrix
    acc = np.empty((0,3)) 
    
    for j in range (0,len(x)):
        ax=0
        ay=0
        az=0
        
        for i in range (0,len(x)):
            if (i!=j):
                dx = x[i]-x[j]
                dy = y[i]-y[j]
                dz = z[i]-z[j]
                norm = np.sqrt(dx**2+dy**2+dz**2)
                ax -= (G*mass[i]*dx)/norm**3
                ay -= (G*mass[i]*dy)/norm**3
                az -= (G*mass[i]*dz)/norm**3
                
        a = np.array([ax,ay,az])
        acc = np.append(acc,[a],axis=0)
        
    return acc


def leapfrog(pos0, vel0, mass, dt):
    
    #-----------------------------------------------------------------
    # INPUT
    # pos0 : initial positions matrix size (N,3). 
    # vel0 : initial velocities matrix size (N,3).
    # mass : masses matrix (1,N).
    # time : time of the simulation ranging from start time (t_start)
    #        to end time (t_end).
    # dt   : timestep.
    #-----------------------------------------------------------------
    # OUTPUT
    # pos_register : dictionnary of all positions matrix size (N,3)
    #                found by leapfrog or Stormer-Verlet integration 
    #                method at each time of the simulation.
    #-----------------------------------------------------------------
    
    global G
    global t_start
    global t_end
    
    pos_register = {}
    
    pos = pos0
    vel = vel0
    acc = NbodyProblem(pos,mass,G)
    
    for t in np.arange(t_start, t_end, dt): # kick-drift-kick version

        vel = vel + acc * dt * 0.5       # kick
        pos = pos + vel * dt             # drift
        
        acc = NbodyProblem(pos,mass,G)   # update acceleration
        vel = vel + acc * dt * 0.5       # kick
        
        pos_register[t] = pos            # save in the register
        
    return pos_register




        
def Coordinates(register, body_id):
    
    #----------------------------------------------------------------
    # INPUT
    # register : Dictionary containing all the positions matrix 
    #            at all time t of the simulation.
    # body_id : identification number of the body, which corresponds 
    #           to the row number in the positions matrix. 
    #----------------------------------------------------------------
    # OUTPUT
    # (Lx,Ly,Lz) : 3 lists containig respectively all x, y and z
    #              values of the desired body (id_body).
    #----------------------------------------------------------------
    
    global t_start
    global t_end
    global dt
    
    Lx = []
    Ly = []
    Lz = []
    for t in np.arange(t_start, t_end, dt):
        matrix = register[t]
        Lx.append(matrix[body_id,0])
        Ly.append(matrix[body_id,1])
        Lz.append(matrix[body_id,2])
        
    return (Lx,Ly,Lz)    
#=======================================================================================================
#                                           Simulation parameters
#=======================================================================================================

N = 5                                                               # number of celestial bodies

t_start = spice.spiceypy.str2et('2022 May 01 00:00:00.0000')        # start time of the simulation

t_end = spice.spiceypy.str2et('2050 May 01 00:00:00.0000')          # end time of simulation

dt = 24*3600                                                           # timestep

G = 6.67430e-11*1e-9                                                    # Newton's Gravitational Constant

num_ts = int(np.ceil(t_end/dt))                                     # number of timesteps

#=======================================================================================================


# EPHEMERIS 2022 May 01 00:00:00.0000

# Sun
sun_id = 0
pos_sun = Ephemeris('SUN')[0]
vel_sun = Ephemeris('SUN')[1]
m_sun = 1.989e30


# Mercury 
mercury_id = 1
pos_mercury = Ephemeris('MERCURY')[0]
vel_mercury = Ephemeris('MERCURY')[1]
m_mercury = 0.330e24


# Mars
mars_id = 2
pos_mars = Ephemeris('MARS BARYCENTER')[0]  
vel_mars = Ephemeris('MARS BARYCENTER')[1]   
m_mars = 0.642e24


# Earth
earth_id = 3
pos_earth = Ephemeris('EARTH')[0]  
vel_earth = Ephemeris('EARTH')[1]  
m_earth = 5.97e24

# Venus
venus_id = 4
pos_venus = Ephemeris('VENUS')[0]    
vel_venus = Ephemeris('VENUS')[1]  
m_venus = 4.87e24


# Initial conditions mass, positions and velocities with ephemeris

M = np.array([m_sun,m_mercury,m_mars,m_earth,m_venus])

P = np.concatenate((pos_sun,pos_mercury,pos_mars,pos_earth,pos_venus),axis=0)

V = np.concatenate((vel_sun,vel_mercury,vel_mars,vel_earth,vel_venus),axis=0)

V -= np.mean(np.dot(M,V),0)/np.mean(M)


# Integration using the leapfrog integration (Störmer-Verlet method)

reg = leapfrog(P,V,M,dt)
print(reg)

# Generating lists containing x,y,z position values of each bodies
Lx_sun = Coordinates(reg,0)[0]
Ly_sun = Coordinates(reg,0)[1]
Lz_sun = Coordinates(reg,0)[2]
print(Lx_sun)

Lx_mercury = Coordinates(reg,1)[0]
Ly_mercury = Coordinates(reg,1)[1]
Lz_mercury = Coordinates(reg,1)[2]

Lx_mars = Coordinates(reg,2)[0]
Ly_mars = Coordinates(reg,2)[1]
Lz_mars = Coordinates(reg,2)[2]

Lx_earth = Coordinates(reg,3)[0]
Ly_earth = Coordinates(reg,3)[1]
Lz_earth = Coordinates(reg,3)[2]

Lx_venus = Coordinates(reg,4)[0]
Ly_venus = Coordinates(reg,4)[1]
Lz_venus = Coordinates(reg,4)[2]


# Ploting the results using Matplotlib

fig = plt.figure("Solar System simulation - Marwa Kadhem")

ax = fig.add_subplot(111, projection='3d')
ax.set_axis_off()
fig.set_facecolor('black')
ax.set_facecolor('black') 

ax.scatter(0,0,0, s=100, c="yellow",label='Sun')
ax.plot(Lx_mercury,Ly_mercury,Lz_mercury, label='Mercury')
ax.plot(Lx_mars,Ly_mars,Lz_mars, label='Mars')
ax.plot(Lx_earth,Ly_earth,Lz_earth, label='Earth')
ax.plot(Lx_venus,Ly_venus,Lz_venus, label='Venus')
ax.legend()

plt.show()

**Tags:** solar-system, orbital-mechanics, python, n-body-simulations
> **Comments:**
> **username_1**: First guess is to make sure all of your units are correct.  This sort of result is common when the force of gravity far exceeds the initial velocities.  
> **username_2**: In addition, numerical stability is much easier if you use AU-solar mass-years units, instead of SI units.  For example, in these units G=4pi.  
> **username_2**: finally, for testing purposes, it would be convenient to break the dependency on spice and  de440 (as these are large file dependencies.)  You are only using them to initialise the system.  So I'd suggest working out the initial values for may 2022 and hard coding them into the program.  When everything is working with the hard-coded values, you can easily switch back the de400 values.  
> **username_3**: Why do you have `G = 6.67430e-11*1e-9`. Where does that extra `*1e-9` come from? Everything else seems to be in mks units but this value. Unrelatedly, I'd suggest you post your code on something like [code review](https://codereview.stackexchange.com/). I think you could really benefit from someone reviewing your code from a coding perspective (aside from the science perspective you're asking for here).  
> **username_4**: Straight lines intersecting the Sun mean that gravity is too strong or speed is too low. A quick way to check is to calculate the *reduced* kinetic and potential energies $v^2/2$ and $-GM_{Sun}/r^2$ and see if the kinetic is about half of the absolute value of the potential energy since a bound orbit will have negative total energy. Next, if the lines radiate *away* from the Sun then you might be missing a minus sign somewhere. Also be extremely careful to never mix kilograms and seconds with kilometers, use only meters (perhaps the 1E-09 zephyr mentions comes from that?)  
> **username_4**: Finally, back up and first see if your method can solve a very familiar one-body problem. In a much shorter, simpler script, just put the Earth at $1.5 \times10^{11} \ \mathbf{\hat{x}}$ meters with a velocity of $2.9780 \times10^4 \ \mathbf{\hat{y}}$ meters per second, and confirm it goes around the $\mathbf{\hat{z}}$ axis  in roughly a circle and takes about a year to do so. All that said, here's a complicated n-body script that matches a JPL ephemeris pretty well https://space.stackexchange.com/a/23409/12102  
> **username_4**: fyi "one-body" means just calculate for the Earth and leave the Sun's central force at the origin. Once you get the integrator working and generating real, physical results, then work your way up slowly; add several bodies without interactions (n x 1-body) then finally add back all the $n(n-1$ two-body interactions i.e. n-body.  
> **username_5**: Since you're doing the calculations yourself, you might want to look at https://astronomy.stackexchange.com/questions/2416/where-can-i-find-a-set-of-data-of-the-initial-conditions-of-our-solar-system/13491#13491 and perhaps https://astronomy.stackexchange.com/questions/13488 (under "If you want to compute star/planet positions yourself, you have several options")  
> **username_0**: Hi everyone, I can't thank you enough for all your suggestions and help. I manage to solve the problem by following the steps you described and by taking into account all your replies.  
> **username_0**: If you want to check the result, it is fully available at https://github.com/MarwaKadhem/Solar-System-3D-simulation  
> **username_0**: And for the demo video: https://m.youtube.com/watch?v=g7ucAwWma2Y  

---

## 0 answers

