2 body n-body program with leapfrog integrator FORTRAN 90

[Joshius43]

So I am having trouble getting this program to work.
I have all of my variables declared and what not, the issue seems to be that when I run the program, the "star" just drifts with an always increasing acceleration continuously away, so the energy of the system is not conserved, as it should be, but rather it is ever increasing.

I have been using print statements to try to identify the problem, but I am unable to see it. Perhaps someone here could help me out?

do n=1, nmax, dt

print *,"-------------------Page Break-------------------"

! Calculate Force
Fr = (G*SM*DMM) / (SQRT(Sx**2+Sy**2)**2)
Frx = cos(Fr)
Fry = sin(Fr)
print *,"Force in x:",Frx,"Force in y:",Fry

!Accelerations
ax = Frx/Sm
ay = Fry/Sm
print *, "Accel in x:",ax," Accel in y:",ay

! Calculate the Velocities in x and y directions with a half timestep
Svx = Svx + 0.5*dt*(ax)
Svy = Svy + 0.5*dt*(ay)
print *, "Speed 1 in x:",Svx," Speed 1 in y:",Svy

! Calculate the New X and Y positions with a full time step
Sx = Sx + dt*Svx
Sy = Sy + dt*Svy
print *,"X position:",Sx," Y Position:", Sy

!Find the new Velocities
Svx = Svx + (dt/2.)*(ax)
Svy = Svy + (dt/2.)*(ax)
print *,"Speed 2 in x:",Svx,"Speed 2 in y:", Svy


!Energy, is it conserved?
E = 0.5*(Sm)*(sqrt(Svx**2+Svy**2))**2
print *, "Energy at ",n," seconds is ", E

if ( n >= nmax) then
print *, "DONE!"
exit
endif

enddo


Thanks in advance

 

[FJacq]

I don't see the potential energy in your energy balance. The kinetic energy increases when the body distance decreases but the potential energy decreases in //. This one is inversely proportional to the distance.

EP = - (G*SM*DMM) / SQRT(Sx**2+Sy**2)

In addition, what means Frx = cos(Fr) ???

Fr is not an angle but a force ! Could you please replace that by something like :

Frx= - Fr * Sx / SQRT(Sx**2+Sy**2)
Fry= - Fr * Sy / SQRT(Sx**2+Sy**2)

Notice the minus sign because the force tends "to reduce" the distance.

A last remark : you might replace (SQRT(Sx**2+Sy**2)**2) by (Sx**2+Sy**2) in the computation of the force

François Jacq

 

[Joshius43]

Thanks FJacq

It turns out that you were right in both cases. I ended up getting it to plot and everything now. As well as implementing an array for my timesteps to allow for far more accuracy. It now simulates the orbit correctly, conserves energy, works. My next step will be moving from this one body orbiting system to a galaxy.

Again, Thanks

Josh Taylor