Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Post Newtonian Gravity Simulation Help

  1. Feb 1, 2015 #1
    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.

    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)

     
     

    Attached Files:

    • gr.png
      gr.png
      File size:
      59.6 KB
      Views:
      78
  2. jcsd
  3. Feb 1, 2015 #2
    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.
     
  4. Feb 2, 2015 #3

    PeterDonis

    User Avatar
    2016 Award

    Staff: Mentor

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Post Newtonian Gravity Simulation Help
Loading...