N-Body Simulation

Setup

Constants

Newton's Version of Kepler's Third Law

Applying to the Earth-Sun System

Celestial Bodies

Sun (☉)

Mercury (☿)

Venus (♀)

Earth (♁)

Mars (♂)

Jupiter (♃)

Saturn (♄)

Uranus (♅)

Neptune (♆)

Pluto (♇)

Indices

Momentum and Energy

Energy

Momentum

Advance Simulation

Advance velocities

Advance positions

N-Body Simulation

The program simulates the gravitational interactions between celestial bodies using numerical methods.

Setup

This section initializes the constants, defines the solar system's parameters, and organizes the celestial bodies into a matrix.

Constants

We start by defining the constants used in the simulation:

days-per-year:=365.24π:=3.141592654solar-mass:=4*π^2Δt:=0.01--

This is equivalent to 1/Δt Hz.

Symbol

Description

Value

days-per-year

Used to normalize planetary data to simulation scale

days-per-year

π

Approximation of pi, critical for calculating orbital mechanics

π

solar-mass

Normalized mass of the Sun, used as a reference unit

solar-mass

Δt

Defines the simulation's time granularity

Δt

The expression solar-mass:=4*π^2 is a result of using normalized astronomical units in conjunction with Kepler's Third Law. This result is derived from the gravitational constant and the mass of the Sun. In normalized units, we set and define the mass of the Sun in such a way that it simplifies calculations in celestial mechanics.

Newton's Version of Kepler's Third Law

Kepler's Third Law, derived from Newtonian mechanics, is:

Solving for the mass of the central body (e.g., the Sun):

Applying to the Earth-Sun System

We can apply this to the Sun-Earth system using astronomical units:

  • AU (astronomical unit)

  • year

  • (solar mass)

Substitute these into the equation:

In astronomical unit systems, it is common to set the gravitational constant for simplicity. Then:

Celestial Bodies

Each celestial body is represented as a vector:

  • Rows 1 to 3: Position (x, y, z).

  • Rows 4 to 6: Velocity (vx, vy, vz).

  • Row 7: Mass (m).

Planets are initialized with a position position, velocity (scaled by days-per-year), and mass (solar-mass fraction). The planets matrix stores the positions, velocities, and masses of all celestial bodies. The source of the data is the JPL HORIZONS celestial mechanics database, and normalized for the solar system.

Sun (☉)

:=[
0
0
0
0
0
0
solar-mass
]
'

Mercury (☿)

Mercury is the closest planet to the Sun and has a highly elliptical orbit.

:=[
-0.38972469318558057 -0.15022403533011131 0.023476059815063703 0.0042923857521344092*days-per-year -0.021425298453769862*days-per-year -0.0023708721826123731*days-per-year 1.660120e-07.*solar-mass
]

Venus (♀)

Venus is the second planet from the Sun and is similar in size and composition to Earth, but has a thick atmosphere and extreme surface conditions.

:=[
-0.71803521454713496 0.041947688018485076 0.041304042456734767 -0.0013897739822126357*days-per-year -0.020119890903497486*days-per-year -0.00030021717325320866*days-per-year 2.447838e-06.*solar-mass
]

Earth (♁)

Earth is the third planet from the Sun and the only known planet to support life. It has a diverse climate and geography.

:=[
-0.17713546150023925 0.96724162210078579 -0.0000039007362829064718 -0.017201146327519165*days-per-year -0.0031864352617521142*days-per-year 0.00000018827814191274*days-per-year 3.003489e-06.*solar-mass
]

Mars (♂)

Mars is the fourth planet from the Sun and is known for its reddish appearance due to iron oxide on its surface. It has a thin atmosphere and polar ice caps.

:=[
1.3907159267594730 -0.013415706106135855 -0.034467796700612273 0.00021345415375975118*days-per-year 0.015123072614264662*days-per-year 0.00030512951698556953*days-per-year 3.227151e-07.*solar-mass
]

Jupiter (♃)

Jupiter is the largest planet in the solar system and is known for its Great Red Spot, a giant storm. It has a thick atmosphere and many moons.

:=[
4.84143144246472090 -1.16032004402742839 -0.103622044471123109 0.00166007664274403694*days-per-year 0.00769901118419740425*days-per-year -0.0000690460016972063023*days-per-year 0.000954791938424326609*solar-mass
]

Saturn (♄)

Saturn is known for its prominent ring system and is the second-largest planet in the solar system. It has a gaseous composition and many moons.

:=[
8.34336671824457987 4.12479856412430479 -0.403523417114321381 -0.00276742510726862411*days-per-year 0.00499852801234917238*days-per-year 0.0000230417297573763929*days-per-year 0.000285885980666130812*solar-mass
]

Uranus (♅)

Uranus is the third-largest planet in the solar system and is known for its blue-green color due to methane in its atmosphere. It has a unique tilt and a faint ring system.

:=[
12.8943695621391310 -15.1111514016986312 -0.223307578892655734 0.00296460137564761618*days-per-year 0.00237847173959480950*days-per-year -0.0000296589568540237556*days-per-year 0.0000436624404335156298*solar-mass
]

Neptune (♆)

Neptune is the eighth and farthest planet from the Sun in the solar system. It is known for its deep blue color and strong winds. Neptune has a dynamic atmosphere and a faint ring system.

:=[
15.3796971148509165 -25.9193146099879641 0.179258772950371181 0.00268067772490389322*days-per-year 0.00162824170038242295*days-per-year -0.000095159225451971587*days-per-year 0.0000515138902046611451*solar-mass
]

Pluto (♇)

Pluto is a dwarf planet located in the Kuiper Belt. It was once considered the ninth planet but was reclassified as a dwarf planet in 2006. Pluto has a highly elliptical orbit and is known for its complex geology and atmosphere.

:=[
-9.87512510193949936 -27.9392880241831424 5.06873275440839440 0.00344178030872987300*days-per-year -0.00152819214839188910*days-per-year -0.00129137458409475460*days-per-year 6.547e-09.*solar-mass
]

Compiled into one matrix:

planets:=[
]
'

Indices

ixes:=combinatorics/n-choose-k(1..=10, 2)ix1:=ixes[1,:]ix2:=ixes[2,:]

The index pairs ix1 and ix2 are used to represent all unique interactions between celestial bodies without repeating any calculations. Each pair of values from ix1 and ix2 identifies two bodies that interact with each other, ensuring that every combination of bodies is considered exactly once.

For example, if one pair specifies that body 1 interacts with body 2, there is no separate calculation for body 2 interacting with body 1, as the interaction is symmetric according to Newton's third law of motion. This means that the gravitational force exerted by body 1 on body 2 is equal and opposite to the force exerted by body 2 on body 1.

Momentum and Energy

This section calculates the momentum and energy of the system, ensuring accurate simulation of physical laws.

~x:=planets[:,[
1
2
3
]
]
~v:=planets[:,[
4
5
6
]
]
m:=planets[:,7]

Squared differences in positions for all body pairs.

Δ:=(x[ix1,:]-x[ix2,:])^2distance:=stats/sum/column(Δ)

Energy

Caculate the system's total energy by summing kinetic and potential energies. Kinetic energy comes from each body's mass and velocity, while potential energy depends on gravitational interactions between body pairs based on their masses and distances. The total energy ensures the simulation conserves energy, reflecting physical accuracy.

kinetic:=0.5*m*stats/sum/column(v^2)potential:=m[ix1]*m[ix2]/distance^0.5ek:=stats/sum/row(kinetic)em:=stats/sum/row(potential)total-energy:=ek-em

Momentum

Calculates the total momentum of the system to ensure the center of mass remains stationary. It computes momentum as the product of each body's mass and velocity, summing across all bodies. An offset is applied to the Sun's velocity to balance the system, maintaining a stable reference frame.

ps:=stats/sum/column(v)*moffset:=ps/-solar-mass

Advance Simulation

Updates the positions and velocities of celestial bodies over time based on gravitational forces.

Advance velocities

Gravitational interactions between body pairs are calculated using their distances and masses. A scaling factor (magnitude) determines how much the velocities change, based on the inverse square of the distance. Velocities are adjusted for each pair, ensuring forces are applied equally and oppositely, maintaining symmetry.

mag:=Δt*distance^-1.5v[ix1,:]-=Δ*m[ix2]*magv[ix2,:]+=Δ*m[ix1]*mag

Advance positions

After updating velocities, positions are recalculated by moving each body according to its velocity and the time step. This step advances the simulation, showing how bodies evolve over time under gravitational influence.

x+=v*Δt