Help Solve FORTRAN Post Newtonian Gravity Simulation

[goalie39]

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.

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)

 

[goalie39]

My apologies, I guess this should have gone in the homework area. If someone could move it, I'd appreciate it. Sorry for the confusion.

 

[PeterDonis]

I guess this should have gone in the homework area.

Not really, since it's not a homework problem, it's just something you're doing in your free time.

However, I'm not sure how many FORTRAN coders there are here, so you might not get much useful input.