Simulating the orbits of the Earth, Moon and Sun together in Fortran

This has been racking my brain for hours now. I will try to explain as much as I can. I'm aware this code is horrendously ungraceful, but I'm new to Fortran and making things more efficient is not easy for me.

I am trying to simulate the celestial motion I've specified in the title. I am going to have (before moving to the co-moving frame) the Sun at the origin of coordinates, the Earth at x = 149597870d+3meters away (the mean distance from the Sun to the Earth I believe), y = 0. For the moon, I'm having it at the distance from the Moon to Earth (384400d+3 meters) plus the distance from the Earth to the Sun, so that we essentially have a situation where all planets are situated on the y = 0 line, with the Earth an astronomical unit away in the x-axis, and the moon a further distance away equal to the distance between the Earth and the moon.

I've attempted to illustrate the situation, as well as the free body diagrams, in the following two images.

From there, I've defined a 12 dimensional array, sol, which holds the x and y positions and velocities of each body such that sol = (x1,y1,x2,y2,x3,y3,dx1/dt,dy1/dt,dx2/dt,dy2/dt,dx3/dt,dy3/dt). I then initialized the array: the sun will have no initial velocities, the Earth will have an initial y velocity equal to its mean orbital velocity around the sun with no initial x velocity, and the moon will have an initial y velocity equal to its mean orbital velocity around the earth and no initial x velocity.

 inits(1:2) = 0d0
inits(3) = distE
inits(4) = 0d0
inits(5) = distE+distM
inits(6:9) = 0d0
inits(10) = vE
inits(11) = 0d0
inits(12) = vM
sol = inits


I then attempt to bring everything to a co-moving frame based off the following equation:

 mf = 0
Mtot = mSun + mEarth + mMoon
mf(1:2) = mSun/(Mtot) * mf(1:2) + mEarth/(Mtot) * sol(3:4) + mMoon/(Mtot) * sol(5:6)
mf(3:4) = mEarth/(Mtot) * mf(3:4) + mSun/(Mtot) * sol(1:2) + mMoon/(Mtot) * sol(5:6)
mf(5:6) = mMoon/(Mtot) * mf(5:6) + mSun/(Mtot) * sol(1:2) + mEarth/(Mtot) * sol(3:4)
mf(7:8) = mSun/(Mtot) * mf(7:8) + mEarth/(Mtot) * sol(9:10) +  mMoon/(Mtot) * sol(11:12)
mf(9:10) = mEarth/(Mtot) * mf(9:10) + mSun/(Mtot) * sol(7:8) + mMoon/(Mtot) * sol(11:12)
mf(11:12) = mMoon/(Mtot) * mf(11:12) + mSun/(Mtot) * sol(7:8) + mEarth/(Mtot) * sol(9:10)
sol = inits - mf


I then use a varied timestep for my calculations based off the equation:

real*8 function fun_tstepq8(arr)
use importants
implicit none
real*8,dimension(12), intent(in) :: arr
real*8::alp,part1,part2,part3,distEtoS,distEtoM,distStoM
real*8 :: distEdottoMdot,distEdottoSdot,distSdottoMdot
alp = 1d-4
distEtoS = SQRT((sol(3)-sol(1))**2 + (sol(4)-sol(2))**2)**3
distEtoM = SQRT((sol(5)-sol(3))**2 + (sol(6)-sol(4))**2)**3
distStoM = SQRT((sol(5)-sol(1))**2 + (sol(6)-sol(2))**2)**3
distEdottoSdot = SQRT((sol(9)-sol(7))**2 + (sol(10)-sol(8))**2)
distEdottoMdot = SQRT((sol(11)-sol(9))**2 + (sol(12)-sol(10))**2)
distSdottoMdot = SQRT((sol(11)-sol(7))**2 + (sol(12)-sol(8))**2)
part1= distEtoS/distEdottoSdot
part2= distEtoM/distEdottoMdot
part3= distStoM/distSdottoMdot
fun_tstepq8 = alp * MIN(part1,part2,part3)
end function


Putting that all together, I use a Forward Euler method and calculate the function f from y0 = y0 + hf using the subroutine:

subroutine q8RHS(sol)
use importants, ONLY: mSun, mEarth,mMoon, F, G
implicit none
real*8 :: distEtoS,distEtoM,distStoM
real*8,dimension(12) :: sol
integer :: i
distEtoS = SQRT((sol(3)-sol(1))**2 + (sol(4)-sol(2))**2)**3
distEtoM = SQRT((sol(5)-sol(3))**2 + (sol(6)-sol(4))**2)**3
distStoM = SQRT((sol(5)-sol(1))**2 + (sol(6)-sol(2))**2)**3
do i = 1,12
if (i < 7) then
F(i) = sol(i+6)
elseif (i < 9) then
F(i) = G * (mEarth * (sol(i-4) - sol(i-6))/distEtoS + mMoon * (sol(i-2) - sol(i-6))/distStoM)

elseif (i < 11) then
F(i) =  G * mSun * (sol(i-6) - sol(i-8))/distEtoS - mMoon * G *(sol(i-4) - sol(i-6))/distEtoM
else
F(i) =  -G * mSun * (sol(i-6) - sol(i-10))/distStoM -G*mEarth * (sol(i-6) - sol(i-8))/distEtoM

endif
enddo
end subroutine


In terms of deriving the q8RHS subroutine, I started with:

And then this is my work to derive mechanics for one of the planets. Note I use the Euler method here, and derive suitable functions to use for f(x,t) = dx/dt form needed for the Euler method. Note I couldn't use the denominator I derived as at some point x2 = x1 in my simulation if I use this, hence the issue I talked about with the timestep getting tiny enough where the simulation grinds to a halt.

Overall, my program looks like this:

module importants
implicit none
integer,parameter::Ndim = 3
integer:: N,N2
real*8, parameter :: G = 6.67408e-11
integer,parameter :: Nbodies = 2
integer::specific_Nbodies = 2*Nbodies
real*8,dimension(12) :: inits,sol,F
real*8 :: d_j = 778.547200d+9
real*8 :: v_j = 13.1d+3
integer :: rank = Nbodies * Ndim * 2
real*8 :: mSun = 1.989d+30
real*8, dimension(12) :: mf
real*8 :: mJup = 1.89819d+27
real*8, parameter :: pi= DACOS(-1.d0)
real*8 :: a,P,FF
real*8 :: pJ = 374080032d0
real*8, dimension(2) :: QEcompare,QLcompare,QPcompare
real*8 :: mEarth = 5.9722d+24
real*8 :: mMoon = 7.34767309d+22
real*8 :: distE = 149597870d+3 ! Dist from sun to earth
real*8 :: distM = 384400d+3 ! Dist from earth to moon
real*8 :: distMS = 152d+9
real*8 :: vE = 29.78d+3
real*8 :: vM = 1.022d+3
end module

program exercise3
use importants
implicit none
print*,'two'
!call solve()
!call q3()
!call q4()
!call q5()
!call q5circ()
!call q6ecc()
!call q6pt2()
!call q7()
call q8()
end program

subroutine q8()
use importants
implicit none
!real*8,external :: Epot,Ekin, L
real*8,external :: fun_tstepq8
real*8 :: tstep,time,tracker,Mtot,t
!real*8,dimension(2) :: Etot,L1,L2
!real*8, dimension(3,2) :: r2
integer :: j,i
print*, 'Starting Q8. -------------------------------------------------------'
inits(1:2) = 0d0
inits(3) = distE
inits(4) = 0d0
inits(5) = distE+distM
inits(6:9) = 0d0
inits(10) = vE
inits(11) = 0d0
inits(12) = vM
sol = inits
mf = 0
Mtot = mSun + mEarth + mMoon
mf(1:2) = mSun/(Mtot) * mf(1:2) + mEarth/(Mtot) * sol(3:4) + mMoon/(Mtot) * sol(5:6)
mf(3:4) = mEarth/(Mtot) * mf(3:4) + mSun/(Mtot) * sol(1:2) + mMoon/(Mtot) * sol(5:6)
mf(5:6) = mMoon/(Mtot) * mf(5:6) + mSun/(Mtot) * sol(1:2) + mEarth/(Mtot) * sol(3:4)
mf(7:8) = mSun/(Mtot) * mf(7:8) + mEarth/(Mtot) * sol(9:10) +  mMoon/(Mtot) * sol(11:12)
mf(9:10) = mEarth/(Mtot) * mf(9:10) + mSun/(Mtot) * sol(7:8) + mMoon/(Mtot) * sol(11:12)
mf(11:12) = mMoon/(Mtot) * mf(11:12) + mSun/(Mtot) * sol(7:8) + mEarth/(Mtot) * sol(9:10)
sol = inits - mf
print*, '---------------------------------------------------------------------'
print*, 'Beginning fixed t timestep. '
t = 0d0
P = 3.154d+7
time=0
tracker=0
print*, 'P understood as :', P
open(70,file='q8.dat')
time=0
tracker=0
write(70,*) sol(1),sol(2),sol(3),sol(4),sol(5),sol(6)
do while(time.le.P)
call q8RHS(sol)
tstep=fun_tstepq8(sol)
time=time+tstep
tracker=tracker+1
print*, tstep
do i=1,12
sol(i) = sol(i) + tstep * F(i)
enddo
write(70,*) sol(1),sol(2),sol(3),sol(4),sol(5),sol(6)

enddo

end subroutine


I know there are a bunch of variables that are not used here, but were used for other subroutines, as this is a bit, ungainly beast of a program I've made. I just know that the previous subroutines and all of their related variables have run fine, and I have no issue with them and they don't interfere with the code I have in mind.

The motions of the bodies look like this:

I'm not sure if the Earth and Moon are behaving properly with regards to each other. Zooming in, the moon does seem to oscillate to the left and right of the Sun's trajectory, so if the z axis was viewable this might make more sense, but it also might not.

The behavior above of the Moon's path moving to the left and right of the Earth's repeats constantly, but I'm not sure if this is still necessarily physically sensible.

• CodeReview is for working code only, to improve clarity, efficiency, and maintainability. It does not appear your code is working yet (but it might be working fine with poor initial conditions, difficult to say), so is off-topic here. Try StackOverflow. Apr 11, 2020 at 17:32
• I haven't study your code in depth, but "the moon will have an initial y velocity equal to its mean orbital velocity around the earth" is definitely wrong. In a fixed reference frame, the moon's initial velocity should be "mean orbital velocity around earth plus earth mean orbital velocity around the sun".
– vnp
Apr 11, 2020 at 19:55
• @AJNeufeld I don't think (yet) that this code doesn't work. The results seem to be consistent with the initial conditions. The fact that the conditions are wrong does not invalidate the code.
– vnp
Apr 11, 2020 at 19:59
• Please can you translate your images to MathJax. Note that we use $ rather than  for inline formula. $f$ ->$f\\$. This increases the accessibility of your post for users with visual impairments. Apr 12, 2020 at 0:32
• Thanks for your efforts, and please keep in mind the code should work before posting the first version next time.
– Mast
Apr 12, 2020 at 14:12

This question is an oldie but goodie. Nice work! It really deserves some love. An edit bumped it recently to the front page, and that got me interested in how I’d refactor the program.

The Modern Conveniences

You declare your variables in the style of Fortran-77, which is fine. Twenty-first-century Fortran has some new syntax that will improve your quality of life by letting you write structured code instead of labels, simplify your array declarations, and (in theory) have full standardized portability for your code. For example, the declaration of the masses of the celestial bodies, currently three mutable variables, might become the array:

real(real64),parameter :: masses(nBodies) = [ 1.989d+30, 5.9722d+24, 7.34767309d+22 ]


Don’t Let Anything Modify Your Data that Doesn’t Need to

The more places in your code can modify a piece of data, the harder it is to figure out where a modification came from. The compiler will also be able to optimize better, especially when compiling a loop that calls a function or subroutine, if it knows nothing could have changed a variable behind its back.

One corollary of this is that you should declare parameter data wherever you can. The code currently does this for some vonstants, such as G, but not others.

Another aspect of this is that it makes sense to factor your program into small functions and subroutines that each declare only the local data they need. There are fewer things to keep track of and fewer variables you could mix up by mistake.

Use Multidimensional Arrays

Currently, you do most of your work on flat arrays that you access by index, with expressions such as:

mf(7:8) = mSun/(Mtot) * mf(7:8) + mEarth/(Mtot) * sol(9:10) +  mMoon/(Mtot) * sol(11:12)


This is confusing. You already declare Ndim and Nbodies as integer parameters. You could instead declare your arrays as, for example

real(real64) :: positions(nBodies,nDim)
real(real64) :: velocities(nBodies,nDim)


If you then remember that each column holds the data for a celestial body, in order of decreasing mass, and each row holds the coordinates for an individual body, it becomes easy to see that velocity(i,1) is the x-coordinate of celestial body i. You can still define parameter constants for the indices of the sun, earth and moon if you like. (But I advise hardcoding those into your program logic; I’ll come back to this later.)

Use Array Operations

These are not only very easy to write and to read, they optimize very well. Let’s suppose we have multidimensional arrays holding the component-wise positions, velocities and accelerations of each vector. A simple subroutine, using only first-year calculus, to update these arrays, would look like this:

subroutine update_system( dt, ps, vs, as )
use importants
implicit none
real(real64) ,intent(in) :: dt
real(real64) :: ps(nBodies,nDim)
real(real64) :: vs(nBodies,nDim)
real(real64),intent(in) :: as(nBodies,nDim)

ps = ps + dt*(vs + dt*0.5*as)
vs = vs + dt*as
end subroutine


On a modern compiler with optimizations enabled (such as ICX), this will auto-vectorize..

Make Your Algorithm More General, if Possible

Currently, you’re hardcoding a set of formulas that assume you are working with only the sun, earth and moon, and flat arrays mf and sol laid out in a particular order. Now, a lot of the time, when I say something like, “You might later want to switch between the two-dimensional and three-dimensional cases based on the inclination of the orbits from the ecliptic, or work with a different number of bodies,” you should say, YAGNI: You Ain’t Gonna Need It. Here, though, I think the code becomes a lot more readable and maintainable if you roll the algorithm up into a loop that works for any number of bodies.

Here’s a very simple function—using only high-school physics—that calculates each component of acceleration for each body from the current positions and masses of the bodies. All inputs and outputs are arrays with dimension (nBodies, nDim). The elementary round-trip polar conversion only supports two dimensions, but it allows an arbitrary number of bodies, without assumptions about their names or ordering.

pure function accelerations(ps) result(as)
use importants
implicit none
real(real64),intent(in) :: ps(nBodies,nDim)
real(real64) :: as(nBodies,nDim)

real(real64) :: rsquared
real(real64) :: angles(nBodies,nBodies)
real(real64) :: forces(nBodies,nBodies)
integer :: i, j
real(real64) :: components(nDim)

! Not needed in theory, but makes the behavior more deterministic.
do i = 1, nBodies
do j = 1, nBodies
angles(i,j) = 0.0
forces(i,j) = 0.0
end do
end do

! For all i > j, set forces(i,j) to the gravitational force exerted on
! body i by body j, and angles(i,j) to the direction of that vector.
do i = 1, nBodies
do j = 1, i-1
rsquared = (ps(j,1)-ps(i,1))**2 + (ps(j,2)-ps(i,2))**2
angles(i,j) = atan2(ps(j,2)-ps(i,2), ps(j,1)-ps(i,1))
forces(i,j) = G*masses(i)*masses(j)/rsquared
end do
end do

do i = 1, nBodies
components = [0.0, 0.0]
do j = 1, i-1
components(1) = components(1) + (forces(i,j) * cos(angles(i,j)))
components(2) = components(2) + (forces(i,j) * sin(angles(i,j)))
end do

! Take advantage of the symmetry of the problem: the force exerted by j
! on i is equal and opposite to the force exerted by i on j.
do j = i+1, nBodies
components(1) = components(1) - (forces(j,i) * cos(angles(j,i)))
components(2) = components(2) - (forces(j,i) * sin(angles(j,i)))
end do

as(i,1) = components(1) / masses(i)
as(i,2) = components(2) / masses(i)
end do
end function


I make no claim that this is a good implementation. The error is much higher than with a Runge–Kutta method like the one you used. It’s also very plausible there’s a bug in there that I haven’t caught. It seems to behave reasonably with real-world initial conditions for about a month.

Document the Code

You currently have short, often-cryptic names, such as sol(9:10), and very few comments. This makes it very hard for someone else (or myself a year or two later when I’ve done that) to figure out what the code is doing or whether you’re actually using the right slices of sol and mf.

It would be better to put at least some of thee explanation of your algorithm that you currently have as scans of your handwritten notes in the post, as comments in the code itself, or at minimum a URL.

A Minor Issue on Types

You declare most of your variables REAL*8. I don’t know of any compiler in actual use that rejects that but accepts the formal standard, but in principle, the language defines the constant real64 in the module ISO_FORTRAN_ENV as the kind of a 64-bit REAL. I therefore write real(real64), but many programmers abbreviate it. Aliasing this to something like real(realKind) instead would make it easier to change the precision.

Putting it All Together

Here’s what I finally came up with.

module importants
use iso_fortran_env, only:real64
implicit none
integer,parameter :: nDim = 2
! All columns of bodies are in the order [Sun, Earth, Moon].
integer,parameter :: nBodies = 3
real(real64),parameter :: masses(nBodies) = [ 1.989d+30, 5.9722d+24, 7.34767309d+22 ]
real(real64),parameter :: G = 6.67408e-11
end module

module interfaces
interface
subroutine update_system( dt, ps, vs, as )
use importants
real(real64),intent(in) :: dt
real(real64) :: ps(nBodies,nDim)
real(real64) :: vs(nBodies,nDim)
real(real64),intent(in) :: as(nBodies,nDim)
end subroutine

pure function accelerations(ps) result(as)
use importants
real(real64),intent(in) :: ps(nBodies,nDim)
real(real64) :: as(nBodies,nDim)
end function
end interface
end module

! Given a time increment and component-wise position, velocity and
! acceleration vectors of the same shape, advances the position and velocity.
subroutine update_system( dt, ps, vs, as )
use importants
implicit none
real(real64) ,intent(in) :: dt
real(real64) :: ps(nBodies,nDim)
real(real64) :: vs(nBodies,nDim)
real(real64),intent(in) :: as(nBodies,nDim)

ps = ps + dt*(vs + dt*0.5*as)
vs = vs + dt*as
end subroutine

! Returns the component-wise acceleration based on the arrays of the coord-
! inates and masses of the bodies.
pure function accelerations(ps) result(as)
use importants
implicit none
real(real64),intent(in) :: ps(nBodies,nDim)
real(real64) :: as(nBodies,nDim)

real(real64) :: rsquared
real(real64) :: angles(nBodies,nBodies)
real(real64) :: forces(nBodies,nBodies)
integer :: i, j
real(real64) :: components(nDim)

! Not needed in theory, but makes the behavior more deterministic.
do i = 1, nBodies
do j = 1, nBodies
angles(i,j) = 0.0
forces(i,j) = 0.0
end do
end do

! For all i > j, set forces(i,j) to the gravitational force exerted on
! body i by body j, and angles(i,j) to the direction of that vector.
do i = 1, nBodies
do j = 1, i-1
rsquared = (ps(j,1)-ps(i,1))**2 + (ps(j,2)-ps(i,2))**2
angles(i,j) = atan2(ps(j,2)-ps(i,2), ps(j,1)-ps(i,1))
forces(i,j) = G*masses(i)*masses(j)/rsquared
end do
end do

do i = 1, nBodies
components = [0.0, 0.0]
do j = 1, i-1
components(1) = components(1) + (forces(i,j) * cos(angles(i,j)))
components(2) = components(2) + (forces(i,j) * sin(angles(i,j)))
end do

! Take advantage of the symmetry of the problem: the force exerted by j
! on i is equal and opposite to the force exerted by i on j.
do j = i+1, nBodies
components(1) = components(1) - (forces(j,i) * cos(angles(j,i)))
components(2) = components(2) - (forces(j,i) * sin(angles(j,i)))
end do

as(i,1) = components(1) / masses(i)
as(i,2) = components(2) / masses(i)
end do
end function

program threebody
use importants
use interfaces
implicit none

real(real64),parameter :: tn = 60.d0*60*24*28
real(real64),parameter :: deltat = 1
real(real64) :: t = 0
real(real64) :: positions(nBodies,nDim)
real(real64) :: velocities(nBodies,nDim)

! Set the initial conditions in column-major order.
positions = reshape(([0.d0, 1.47095d+11, 1.47095d+11 + 4.055d+8, 0.d0, 0.d0, 0.d0]), shape(positions))
velocities = reshape(([0.d0, 0.d0, 0.d0, 0.d0, 2.978d+3, 2.978d+3 + 970.d0]), shape(velocities))

print *, positions(1:1, :), positions(2:2, :), positions(3:3, :)
print *, sqrt((positions(2,1)-positions(1,1))**2 + (positions(2,2)-positions(1,2))**2), sqrt((positions(3,1)-positions(2,1))**2 + (positions(3,2)-positions(2,2))**2)

do while (t < tn)
call update_system( deltat, positions, velocities, accelerations(positions) )
!print *, positions(1:1, :), positions(2:2, :), positions(3:3, :)
t = t + deltat
end do

print *, positions(1:1, :), positions(2:2, :), positions(3:3, :)
print *, sqrt((positions(2,1)-positions(1,1))**2 + (positions(2,2)-positions(1,2))**2), sqrt((positions(3,1)-positions(2,1))**2 + (positions(3,2)-positions(2,2))**2)
end program


I’m certain you can improve on the algorithm. The general structure here, and in particular using arrays with the same intuitive shape, and passing them to and returning them from functions, is a strategy you and others can use.

• Yay, bravo! Generalizing beyond 3 bodies is helpful. I especially appreciate improvements around "cryptic names, such as sol(9:10)".
– J_H
Jan 29 at 2:39
program celestial_motion
implicit none

real :: xe, ye, xm, ym, xs, ys, dt, t
real :: G, Me, Ms, mm

! Constants
G = 6.67430d-11 ! gravitational constant
Me = 5.9722d24 ! mass of the Earth
Ms = 1.9885d30 ! mass of the Sun
mm = 7.342d22 ! mass of the Moon

! Initial conditions
xe = 149597870e3 ! initial x-coordinate of Earth
ye = 0.0 ! initial y-coordinate of Earth
xm = xe + 384400e3 ! initial x-coordinate of Moon
ym = 0.0 ! initial y-coordinate of Moon
xs = 0.0 ! initial x-coordinate of Sun
ys = 0.0 ! initial y-coordinate of Sun

! Time step
dt = 60.0 ! 60 seconds
t = 0.0 ! initial time

! Simulation loop
do while (t < 1d7)
! Calculate gravitational force between Earth and Sun
xe = xe + dt * (G * Ms * (xs - xe) / ((xs - xe)**2 + (ys - ye)**2)**1.5)
ye = ye + dt * (G * Ms * (ys - ye) / ((xs - xe)**2 + (ys - ye)**2)**1.5)

! Calculate gravitational force between Earth and Moon
xe = xe + dt * (G * mm * (xm - xe) / ((xm - xe)**2 + (ym - ye)**2)**1.5)
ye = ye + dt * (G * mm * (ym - ye) / ((xm - xe)**2 + (ym - ye)**2)**1.5)

! Calculate gravitational force between Moon and Sun
xm = xm + dt * (G * Ms * (xs - xm) / ((xs - xm)**2 + (ys - ym)**2)**1.5)
ym = ym + dt * (G * Ms * (ys - ym) / ((xs - xm)**2 + (ys - ym)**2)**1.5)

! Update time
t = t + dt
end do
end program celestial_motion


This code implements a simple numerical simulation of the motion of the Earth, Moon, and Sun using the gravitational force between the celestial bodies. The simulation calculates the gravitational force on each body using the Newtonian formula for gravity, and updates the position of each body at each time step. You can adjust the time step (dt) and the length of the simulation (t) to see how the orbits change over time.

• Do you think you could say something about how the OP could improve their code? Answers that only present a different solution, without a review of the original code, are discouraged here. Jan 30 at 0:51
• Also, this code doesn’t appear to be correct. At a glance. It neither calculates the pull of the Earth on the Moon, nor any of the other celestial bodies on the Sun. Neither could it be extended to N bodies or D dimensions without writing out all DN(N-1) cases. making it very difficult to re-use. Jan 30 at 9:52
• Or if these are meant to be in a frame of reference where omitting half the checks works out, the derivation should be explained. Jan 30 at 10:22