# 2-body problem with Fortran 95

I am trying to implement the N-body problem using Fortran 95. Before generalizing it to N bodies, I am first trying to calculate the motion of a single orbiting body (body 2 in my code below) around a fixed one (body 1 in the code below). The problem is the output is not correct, but I can't find the error with my equations! This is what I am trying to implement:

$V_{2,x} \leftarrow V_{2,x} - \frac{ G \times M_1 \times (x_2-x_1)}{r^3} \times dt$

$V_{2,y} \leftarrow V_{2,y} - \frac{ G \times M_1 \times (y_2-y_1)}{r^3} \times dt$

$x_2 \leftarrow x_2 + V_{2,x} \times dt$

$y_2 \leftarrow y_2 + V_{2,y} \times dt$

This is my code:

program gravitation
implicit none

integer :: k, total
real, dimension(2) :: M, x, y, Vx, Vy
real, parameter :: G = 6.6738E-11, pi = 3.1415926
real :: time, dt, r

dt = 1000. ! dt = 1000 seconds
total = 5000 ! total number of points to be output

M(1)=1.9891E+09 ! Sun's mass
x(1)=0.
y(1)=0.
Vx(1)=0.
Vy(1)=0.

M(2)=5.9736E+03 ! Earth's mass
x(2)=1.4960E+11 ! Earth's average orbit as initial point
y(2)=0.
Vx(2)=0.
Vy(2)=2.9786E+04 ! Earth's average orbital speed

open(unit=101,file="a.txt")
write(101,*) "time, x, y"
write(101,"(E10.4,2(A1,1x,E10.4))") 0., ",", x(2), ",", y(2)

do k=1,total
r = sqrt( (x(2)-x(1))**2 + (y(2)-y(1))**2 )
Vx(2) = Vx(2) - G * M(1) * dt * (x(2)-x(1)) / r**3 ! i.e., Vx2 = Vx2 - G.M.dt.(x2-x1)/r^3
Vy(2) = Vy(2) - G * M(1) * dt * (y(2)-y(1)) / r**3
x(2) = x(2) + Vx(2) * dt
y(2) = y(2) + Vy(2) * dt
write(101,"(E10.4,2(A1,1x,E10.4))") k*dt, ",", x(2), ",", y(2)
enddo

end program gravitation


The output file shows that the x position of the second body does not change at all (it is stuck at the initial value of 0.1496E+12), while the y position changes.

Does anyone see what is the problem here? Thanks a lot.

• My suspicion is that your units need work. I haven't checked closely, but I have coded these problems and getting those correct is the first issue to deal with. Try printing your computed acceleration to make sure it is reasonable. Feb 16, 2014 at 8:21
• The usual question applies: What have you done to debug your problem? Which line do you suspect is wrong? Have you run your program in a debugger to inspect why a variable does not change when you expect it to change? Etc. Feb 16, 2014 at 18:25
• @TomDickens I don't think the problem is with my units. I double checked them, and they seem to be all correct (standard SI units). A dt of 1000s seem to be fair as well, and I actually played around with this value but nothing really helps the output... Feb 16, 2014 at 22:35
• I have written working codes like this, and it is standard to use units such that $GM_{sun} = 1$, then the unit of time is the year and the unit of distance is the AU (astronomical unit = mean distance from Earth to Sun). Then you get numbers that are of a reasonable size to work with - for example, JPL specifies state vectors for solar system objects in terms of AU and AU/year. These numbers are of order 1-10 or so... Feb 16, 2014 at 23:23
• @gilberto.agostinho.f: you should write up the fix as an answer and accept it, it'll be easier for users in posterity to see it and learn from the example. Feb 18, 2014 at 23:32

The only problem with the code above was with the values of the masses. By simply correcting the masses of both astronomical bodies to their proper value, the simulation runs pretty well. Below, I added the corrected code (still using all units in SI, i.e., kg, m, seconds); but as one of the users previously commented (TomDickens), probably the best approach for this kind of program would be to use astronomical units, such as AU for distance, year for time and $G \times M_{sun}=1$ concerning masses.

program gravitation
implicit none

integer :: k, total
real(8), dimension(2) :: M, x, y, Vx, Vy
real(8), parameter :: G = 6.6738D-11
real(8) :: time, dt, r

dt = 3600. ! dt = 3600 = 1h
total = 10000 ! total number of points to be output

M(1)=1.9891D+30 ! Sun's mass
x(1)=0.
y(1)=0.
Vx(1)=0.
Vy(1)=0.

M(2)=5.9736D+24 ! Earth's mass
x(2)=1.4960D+11 ! Earth's average orbit as initial point
y(2)=0.
Vx(2)=0.
Vy(2)=2.9786D+04 ! Earth's average orbital speed

open(unit=101,file="2-body_output.txt")
write(101,*) "time, x, y" ! CSV file header
write(101,*) 0., ",", x(2), ",", y(2) ! initial position

do k=1,total
r = sqrt( ( x(2) - x(1) )**2 + ( y(2) - y(1) )**2 ) ! distance between object 1 and 2
Vx(2) = Vx(2) - G * M(1) * dt * ( x(2) - x(1) ) / r**3 ! i.e., Vx2 = Vx2 - G.M.dt.(x2-x1)/r^3
Vy(2) = Vy(2) - G * M(1) * dt * ( y(2) - y(1) ) / r**3
x(2) = x(2) + Vx(2) * dt ! calculating new x position
y(2) = y(2) + Vy(2) * dt
time = k * dt
write(101,*) time, ",", x(2), ",", y(2)
enddo

end program gravitation