Introduction

I remember, when I was an undergraduate student at the University of Texas Rio Grande Valley, I was working with an astrophysics professor, Dr. Pereyra. After teaching me how to use the IDE VIM and the foundation of the programming language Fortran 90 and after seeing my fascination with the simulation of celestial objects, we decided to dive head first at creating an N-Body Simulation program in Fortran 90. We first created code which simulated physics on a single point mass that we set to the sun, then slowly added granularity by generalizing the code. After, making the code general so that point masses could be added dynamically, without compiling, It was possible to take in the coordinate information from the Jet Propulsion Laboratory's Horizons System and simulate the solar system.

The Physics

Newtonian physics was used to simulate the solar system, meaning that the implications of General and Special relativity are not simulated. That being said, by Newtonian physics, I mean using Newton's laws of motion.
The first law states that an object will remain at rest or in uniform motion in a straight line, unless affected by a force. The second law explains that the sum of the forces is equal to the mass of the object multiplied by the acceleration of the object, $\vec{F}=m\vec{a}$. Finally, the third and last law states that when a force is exerted on an object an equal and opposite equal force is exerted. We will use Newton's law of universal gravitation, $\vec{F}_{21}=-G \frac{m_1 m_2}{|r_{12}|^2} \hat{r}$ And by using the second law, $$\vec{g}(\vec{r})=-G\frac{m_1}{|\vec{r}|^2}\hat{r}$$
Now, when we say that the acceleration is constant we mean that the acceleration is equal to some number and doesn't change, $a=c$ and $da=0$. We know from the definition of acceleration that it basically is the change of velocity, $a=\frac{dv}{dt}$. To get the velocity from the acceleration integration is necessary, $$adt=dv\rightarrow a\int_{t_i}^{t_f}dt=\int_{v_i}^{v_f}dv\rightarrow a(t_f-t_i)=v_f-v_i$$ And, setting the change in time constant $t_f-t_i=t$ $$at+v_i=v_f$$ We integrate again to obtain the change in position $$\frac{1}{2}at^2+v_it+x_i=x_f$$ The algorithm then consists of obtaining the initial positions and velocities from a database, calculating the distances between an $i$th object and all the other objects and using those distances to extrapolate the forces and gravitational acceleration, calculating the new velocities and distances, and, finally, repeating as many times as needed.

Code

Fortran 90's structure generally consists of the main "program", subroutines, functions, types, and modules. We will mostly use subroutines and modules. Modules to keep constant global variables and subroutines to easily organize the code. Subroutines allow us to take apart a problem into smaller manageable pieces. Modules is like a package of variables, functions, and subroutines.
Note: I wrote this code when I was first learning how to code Fortran 90, so there are a lot of "slow code". Like how there is no need to calculate the force but you could calculate the acceleration directly.

Modules

We create three different modules and put them in a folder called "modules". The first, "mod_numeric.f90", will contain the kind for integer and real variables. For instance, int_kind will be 4 if you want the integers of your program to be 32 bit and 8 for 64 bit. We will also have the maximum possible number in a kind 8 real or, in other words, a double.

module mod_numeric
implicit none
public
    integer(kind=4), parameter :: int_kind  = 4 ! kind for Integer vars
    integer(kind=4), parameter :: real_kind = 8 ! kind for Real vars
    real(kind=8),    parameter :: max_general = 1.7D308
    real(kind=8),    parameter :: min_general = -1.7D308
end module mod_numeric

The second, "mod_math.f90", will contain a couple of constants

module mod_math
use mod_numeric, only: int_kind, & ! kind for Integer vars
                       real_kind   ! kind for Real var
implicit none
public
    real(kind=real_kind), parameter :: zero = 0D0 ! Number Zero
    real(kind=real_kind), parameter :: one  = 1D0 ! Number One
    real(kind=real_kind), parameter :: two  = 2D0 ! Number Two
    real(kind=real_kind), parameter :: onehundred  = 1D2 ! Number 100
    real(kind=real_kind), parameter :: onethousand = 1D3 ! Number 1000
end module mod_math

The third and last, mod_astronomy.f90, will contain astronomical constants like the gravitational constant, the mass of the sun, etc...

module mod_astronomy
use mod_numeric, only: int_kind, & ! kind for Integer vars
                       real_kind   ! kind for Real vars
implicit none
public
    real(kind=real_kind), parameter :: G = 6.671D-11
    real(kind=real_kind), parameter :: msun = 1.9891D30
    real(kind=real_kind), parameter :: mearth = 5.9742D24
    real(kind=real_kind), parameter :: AU = 1.4959787066D11
    real(kind=real_kind), parameter :: rsun = 6.955D08
    real(kind=real_kind), parameter :: km = 1D3
    real(kind=real_kind), parameter :: Megam = 1D6
    real(kind=real_kind), parameter :: KPS  = 1D3
    real(kind=real_kind), parameter :: CMPS = 1D-2
    real(kind=real_kind), parameter :: day = 86400D0
    real(kind=real_kind), parameter :: hour = 3600d0
    real(kind=real_kind), parameter :: minute = 60d0
end module mod_astronomy

Note: the explanation for the use of the variables are in the source code.

Subroutines

calcforces

The most important subroutine that will be created is the one which will calculates the forces between the point masses. The subroutine which will be named "calcforces" will take in the multidimensional arrays which contain the masses and positions of the point masses, it will also take in the number of point masses. It will have an output multidimensional array which contain the calculated forces. Also, we will have an option to have a fixed point mass at origin of euclidean space along with it's mass also as an input, for testing purposes.

Inputs:

ef_type : 0 - no external force, 1 - point mass at origin
ef_mass : Mass of Point Mass if ef_type = 1
N : Number of point masses
mass : Mass of point masses (will be a 1D array of size N)
pos : Positions of point masses in meters (will be 2D array of size $3\times N$)

Outputs:

forces : Calculated Forces (will be 3D array of size N)

We will have 4 real local arrays with various functions, 2 integer index variables, and 3 real regular local variables. We will also need to import the gravitational constant and the mass of the sun from mod_astronomy and the number zero and one from mod_math.

use mod_math, only: zero, &
                    one
use mod_astronomy, only: G, &
                         msun

The 4 local arrays will have size 3.

real(kind=real_kind) :: position1(1:3), &
                        position2(1:3), &
                        relpos(1:3),    &
                        univec(1:3)

We will name the 2 index variable $i$ and $j$.

integer(kind=int_kind) :: i, j

The three regular variables will be named as follows.

real(kind=real_kind) :: magpos, magpos2, magfor

For the implementation, Firstly, we need to calculate the external forces. We will use the intrinsic function dot_product to calculate the dot product of the vectors.

if (ef_type == 0) then
  forces = zero
else if (ef_type == 1) then
  do i = 1, N
    realpos = pos(1:3,i)
    magpos2 = dot_product(relpos, relpos)
    magpos = sqrt(magpos2)
    univec = relpos / magpos
    forces(1:3,i) = - G * ef_mass * mass(i) &
                        / magpos2 * univec
  end do
end if

Now, to calculate the forces between point masses

do i = 1, N
  position1(1:3) = pos(1:3,i)
  do j = 1, N
    if (j == i) cycle
    position2(1:3) = pos(1:3,j)
    relpos = position1 - position2
    magpos2 = dot_product(relpos, relpos)
    magpos  = sqrt(magpos2)
    forces(1:3,i) = forces(1:3, i)          &
                    - G * mass(i) * mass(j) &
                    / magpos2 * univec
  end do
end do

Finally, return and end ...

  return
end subroutine calcforces

Evolve

Let's call the subroutine that progresses time a step, evolve. It will use the subroutine the calcforce subroutine then use that information to calculate the acceleration, and the new positions and velocities.

Inputs:

ef_type : 0 - no external force, 1 - point mass at origin
ef_mass : Mass of Point Mass if ef_type = 1
deltat : Time step in seconds
N : Number of point masses
mass : Masses of the point masses (array size N)
Input & Outputs
pos : Position of point masses in meters (array size $3\times N$)
vel : Velocity of point masses in meters per second (array size $3\times N$)

It will have local arrays size $3\times N$ for both the forces and acceleration and a local variable $i$ for iteration. Firstly, we call the calcforces subroutine.

call calcforces(ef_type, ef_mass, N, mass, pos, forces)

Then calculate the acceleration

do i = 1, N
  accel(1:3,i) = forces(1:3,i) / mass(i)
end do

Lastly calculate the new positions and the new velocities

pos = pos + deltat * vel + (deltat * deltat) * accel / two
vel = vel + deltat * accel

The compiler will take care of iterating through the arrays element wise.

fimport

We will need to import the data from a text file which contains the x, y, and z elements of the position and velocities of the point masses. The format for the text will be a line base one where the first line is the number of point masses, the next is the the index of the object, 1, then the mass, then the positional x, y, and z, separated by commas, and finally the x, y, and z of the velocity of the object. It then iterates through the rest of the objects. This subroutine will take in the number by which the file was opened, the number of particles. It will output the mass, position, and velocity of the objects.

Inputs:

nunit_particle : the number by which the file was opened
N : the number of point masses
Outputs:
mass : the mass array (size $N$)
pos : the position array (size $3\times N$)
vel : the velocity array (size $3\times N$)

We will have three local integers, one, $i$, will iterate through the lines, i_Object will iterate through the objects, and, Nin, will import the number of objects from the text file. You can see the source code for this subroutine in the link below.

fexport

fexport will output the results in a similar fashion to fimport. It will have a format where the current time will be printed on a text file, then it will iterate through the masses and positions to print them onto the text file. Again, see source code for details.

Putting it all together

The main program

The main program we will call runNBody where we will import a couple of use statements.

program runNBody
use mod_numeric, only: int_kind, & ! kind for Integer vars
                       real_kind   ! kind for Real vars
use mod_math, only: zero, &        ! Number 0
                    onethousand
use mod_astronomy, only: day, &    ! A day in Seconds
                         AU,  &    ! A Astronomical Unit in meters
                         KPS, &    ! A km/s in m/s
                         CMPS, &
                         rsun

Declare the subroutines used as external

external fimport, evolve

Declare the parameters for the numbers used by the Fortran open function that is used to write to text files.

integer(kind=int_kind), parameter :: nunit_indata = 1,   &
                                     nunit_particle = 2, &
                                     nunit_history = 3,  &
                                     nunit_outdata = 4

Declare the maximum number of particle allowed

integer(kind=int_kind), parameter :: Nmax = 999

Declare the variables that will be used

integer(kind=int_kind) :: ef_type
real(kind=real_kind)   :: ef_mass, &
                          deltat
integer(kind=int_kind) :: nTimeSteps, &
                          N,          &
                          file_steps, &
                          scrn_steps

Declare the arrays that will be used

real(kind=real_kind) :: mass(1:N)
real(kind=real_kind) :: pos(1:3,1:Nmax), &
                        vel(1:3,1:Nmax)
real(kind=real_kind)   :: forces(1:3,1:Nmax), &
                          accel(1:3,1:Nmax)

Declare the integers that will be used to iterate

real(kind=real_kind)   :: time
integer(kind=int_kind) :: i, j, k

We will initialize the time to zero

time = zero

and start by reading the text file 'indata.txt' to get the number of particles, the number time steps to take, the size of the time steps, the number of time steps after which to call fexport, the number of time steps after which to output information to the screen, whether to use a point mass at origin, and finally the mass of that point mass.

open( unit = nunit_indata, file = 'indata.t', &
      form='formatted', status='old'
read(nunit_indata,*) N
read(nunit_indata,*) nTimeSteps
read(nunit_indata,*) deltat
read(nunit_indata,*) file_steps
read(nunit_indata,*) scrn_steps
read(nunit_indata,*) ef_type
read(nunit_indata,*)
read(nunit_indata,*) ef_mass
close(nunit_indata)

I output this information for debug purposes in histdata.txt

open( unit = nunit_history, file = 'histdata.txt', &
      form='formatted', status='unknown')
write(nunit_history,*) 'N               = ', N
write(nunit_history,*) 'nTimeSteps      = ', nTimeSteps
write(nunit_history,*) 'deltat          = ', deltat
write(nunit_history,*) 'file_steps      = ', file_steps
write(nunit_history,*) 'screen_steps    = ', scrn_steps
write(nunit_history,*) 'ext. force type = ', ef_type
write(nunit_history,*) 'ext. mass       = ', ef_mass
write(nunit_history,*)
close(nunit_history)

We will now call fimport to import the data and open the file that fexport will use and use it to write the starting data.

call fimport(nunit_particle, N, mass, pos, vel)
open(unit = nunit_outdata, file = 'outdata.txt', &
         form = 'formatted', status = 'unknown')
write(nunit_outdata,*) 'n:(',n,')'
call fexport(nunit_outdata, time, N, mass, pos, vel)

Now for the loop that will iterate through the time steps needed.

do i = 1, nTimeSteps
  call evolve(ef_type, ef_mass, deltat, N, mass, pos, vel, forces, accel
  time = time + deltat
! Prints out status to file every file_steps
  if (i/file_steps*file_steps == i) then
    call fexport(nunit_outdata, time, N, mass, pos, vel)
  end if
! Prints out status to screen every scrn_steps
  if (i/scrn_steps*scrn_steps == i) then
    write(*,"('time=  ', ES21.13, ' Seconds')") time
  end if
end do

Now just close the file and end the program

  close(nunit_outdata)
  stop
end program runNBody

Source Code

Running it

The data file after running is in "FortranNBody/compile/solar system/outdata.txt" in the source code.
In the next part, we will visualize and compute some pretty cool stuff from that data.

Simulating the Solar System : Part 1

Armando Herrera

Armando Herrera

Introduction

The Physics

Code

Modules

Subroutines

calcforces

Evolve

fimport

fexport

Putting it all together

The main program

Running it

Imaging Atmospheric Cherenkov Telescopes/Technique (IACT)

Finding the Higgs boson w/ Shallow Neural Networks

A really brief description of Gravity

C++, keeping Coinbase Pro's level 2 order book in memory

C++, Interacting with Coinbase Pro's Rest API