- #1

- 2

- 0

Hi all:

In my free time, I've been playing with creating a code to help do toy simulations of gravity in my personal code of choice, FORTRAN. The first step was to get Newtonian gravity up and running and this has been pretty successful so far and I've been able to implement about 4 different integrators to play and test with. My next step is to try to get General Relativity going via the Post Newtonian Expansion.

I'm all self taught in this area, so I've been leaning on what I could find by scouring the internet. I found, what seemed to be, a fairly straight forward way to implement this in my code (see the attached image), however, somewhere I've gone wrong and I'm about 2 orders of magnitude off in my calculation. The source of this is here, a website that seems to be dead these days.

I've been staring at this quite a bit and I'm thinking that there has to be something simple I'm just not catching. My hunch is that someone who knows what they're doing might be able to look at it and understand my error fairly quickly. If anyone has the time to read through the attached code, I'd truly appreciate the assistance.

PS - I realize that things are definitely not optimized yet. I'd like to get it working and then focus on making it fast.

Thanks everyone.

In my free time, I've been playing with creating a code to help do toy simulations of gravity in my personal code of choice, FORTRAN. The first step was to get Newtonian gravity up and running and this has been pretty successful so far and I've been able to implement about 4 different integrators to play and test with. My next step is to try to get General Relativity going via the Post Newtonian Expansion.

I'm all self taught in this area, so I've been leaning on what I could find by scouring the internet. I found, what seemed to be, a fairly straight forward way to implement this in my code (see the attached image), however, somewhere I've gone wrong and I'm about 2 orders of magnitude off in my calculation. The source of this is here, a website that seems to be dead these days.

I've been staring at this quite a bit and I'm thinking that there has to be something simple I'm just not catching. My hunch is that someone who knows what they're doing might be able to look at it and understand my error fairly quickly. If anyone has the time to read through the attached code, I'd truly appreciate the assistance.

PS - I realize that things are definitely not optimized yet. I'd like to get it working and then focus on making it fast.

Thanks everyone.

Code:

```
DYDT%DDY = 0D0
DYDT%DDZ = 0D0
C1(:) = 0D0
C2(:) = 0D0
C3(:) = 0D0
A10(:) = 0D0
A1 = 0D0
X_I = BODIES(BODY)%MYPOSITION(1)
Y_I = BODIES(BODY)%MYPOSITION(2)
Z_I = BODIES(BODY)%MYPOSITION(3)
VX_I = BODIES(BODY)%VELOCITY(1)
VY_I = BODIES(BODY)%VELOCITY(2)
VZ_I = BODIES(BODY)%VELOCITY(3)
AX_I = BODIES(BODY)%ACCELERATION(1)
AY_I = BODIES(BODY)%ACCELERATION(2)
AZ_I = BODIES(BODY)%ACCELERATION(3)
M = BODIES(BODY)%MASS
BETA = 1D0
GAMMA = 1D0
A2 = 0D0
DO K = 1,NBODIES
IF(K.EQ.BODY)CYCLE
DX2 = BODIES(K)%MYPOSITION(1) - X_I
DY2 = BODIES(K)%MYPOSITION(2) - Y_I
DZ2 = BODIES(K)%MYPOSITION(3) - Z_I
A2 = A2 + (G_C*BODIES(K)%MASS) / ((DX2*DX2+DY2*DY2+DZ2*DZ2)**(0.5D0))
ENDDO
A2 = 2D0*(BETA+GAMMA)*A2
DO J = 1,NBODIES
IF(J.EQ.BODY)CYCLE
X_J = BODIES(J)%MYPOSITION(1)
Y_J = BODIES(J)%MYPOSITION(2)
Z_J = BODIES(J)%MYPOSITION(3)
VX_J = BODIES(J)%VELOCITY(1)
VY_J = BODIES(J)%VELOCITY(2)
VZ_J = BODIES(J)%VELOCITY(3)
AX_J = BODIES(J)%ACCELERATION(1)
AY_J = BODIES(J)%ACCELERATION(2)
AZ_J = BODIES(J)%ACCELERATION(3)
M_J = BODIES(J)%MASS
!...Distance between bodies
DX = X_J - X_I
DY = Y_J - Y_I
DZ = Z_J - Z_I
A1 = A1 + (G_C*M_J) / ((DX*DX+DY*DY+DZ*DZ)**(1.5D0))
A3 = 0D0
DO K = 1,NBODIES
IF(K.EQ.J)CYCLE
DX2 = BODIES(K)%MYPOSITION(1) - X_J
DY2 = BODIES(K)%MYPOSITION(2) - Y_J
DZ2 = BODIES(K)%MYPOSITION(3) - Z_J
A3 = A3 + (G_C*BODIES(K)%MASS) / ((DX2*DX2+DY2*DY2+DZ2*DZ2)**(0.5D0))
ENDDO
A3 = (2D0*BETA-1D0)*A3
!...Sum of square velocities
A5 = GAMMA*(VX_I*VX_I+VY_I*VY_I+VZ_I*VZ_I)
!...Sum of square velocities - Sum of product of velocities
A6 = (1D0+GAMMA)*((VX_J*VX_J+VY_J*VY_J+VZ_J*VZ_J) - &
2D0*(VX_I*VX_J-VY_I*VY_J-VZ_I*VZ_J))
!...DX*Velocity Term / Square distance
A7 = 1.5D0*((-DX*VX_J-DY*VY_J-DZ*VX_J)**2D0 &
/ (DX*DX+DY*DY+DZ*DZ))
!...DX*Acceleration Term * DX
A8 = 0.5D0*(DX*AX_J+DY*AY_J+DZ*AZ_J)
!...DX * GAMMA * Velocity
A9 = (DX*(2D0+2D0*GAMMA)*VX_I - (1D0+2D0*GAMMA)*VX_J) + &
(DY*(2D0+2D0*GAMMA)*VY_I - (1D0+2D0*GAMMA)*VY_J) + &
(DZ*(2D0+2D0*GAMMA)*VZ_I - (1D0+2D0*GAMMA)*VZ_J)
A10(1) = A10(1) + (G_C*M_J*AX_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))
A10(2) = A10(2) + (G_C*M_J*AY_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))
A10(3) = A10(3) + (G_C*M_J*AZ_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))
!...Combination of terms
C1(1) = C1(1) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DX ) + A9 * (VX_I - VX_J)
C1(2) = C1(2) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DY ) + A9 * (VY_I - VY_J)
C1(3) = C1(3) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DZ ) + A9 * (VZ_I - VZ_J)
ENDDO
DYDT%DDX = ((A1*C1(1)) + 0.5D0*(3D0+4D0*GAMMA)*A10(1))/(C*C)
DYDT%DDY = ((A1*C1(2)) + 0.5D0*(3D0+4D0*GAMMA)*A10(2))/(C*C)
DYDT%DDZ = ((A1*C1(3)) + 0.5D0*(3D0+4D0*GAMMA)*A10(3))/(C*C)
```