https://github.com/florentrenaud/nbody6tt
Raw File
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
fdisk.f
      SUBROUTINE FDISK(XI,XIDOT,FM,FD)
*
*
*       Miyamoto disk force.
*       --------------------
*
      IMPLICIT REAL*8  (A-H,O-Z)
      COMMON/GALAXY/ GMG,RG(3),VG(3),FG(3),FGD(3),TG,
     &               OMEGA,DISK,A,B,V02,RL2,GMB,AR,GAM,ZDUM(7)
      REAL*8  XI(3),XIDOT(3),FM(3),FD(3)
*
*
*       Obtain force & derivative for global variables XI & XIDOT.
      R2 = XI(1)**2 + XI(2)**2
      BZ = SQRT(B**2 + XI(3)**2)
      AZ = SQRT(R2 + (A + BZ)**2)
      AZ3 = DISK/AZ**3
*       Note missing square root sign in (b^2 + z^2) of Book eq. (8.52).
      AZDOT = XI(1)*XIDOT(1) + XI(2)*XIDOT(2) +
     &                         (A + BZ)*XI(3)*XIDOT(3)/BZ
      FM(1) = -AZ3*XI(1)
      FM(2) = -AZ3*XI(2)
      FM(3) = -AZ3*XI(3)*(A + BZ)/BZ
*     RDOT = (XI(1)*XIDOT(1) + XI(2)*XIDOT(2))/SQRT(R2)
      FD(1) = -AZ3*(XIDOT(1) - 3.0*AZDOT*XI(1)/AZ**2)
      FD(2) = -AZ3*(XIDOT(2) - 3.0*AZDOT*XI(2)/AZ**2)
*       Note wrong sign in first term of eq. (8.52) (see printer errata).
      Y1 = 3.0*(A + BZ)*XI(3)*AZDOT/(AZ**2*BZ)
*       Note printer mistake of dZ/dt in denominator (see errata).
      Y2 = (A*B**2 + BZ**3)*XIDOT(3)/BZ**3
      FD(3) = AZ3*(Y1 - Y2)
*       Refer to Miyamoto & Nagai (PASJ 27, 533) and Book eq. (8.52).
*
      RETURN
*
      END
back to top