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

Question

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.

Answer 1

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.

Answer 2

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.