Hi all:(adsbygoogle = window.adsbygoogle || []).push({});

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 (Text):

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)

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Post Newtonian Gravity Simulation Help

**Physics Forums | Science Articles, Homework Help, Discussion**