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
fcorr.f
      SUBROUTINE FCORR(I,DM,KW)
*
*
*       Total force corrections due to masss loss.
*       ------------------------------------------
*
      INCLUDE 'common6.h'
      REAL*8  A(9)
      LOGICAL IKICK
*
*
*       Save the velocity components and square velocity.
      VI2 = 0.0
      DO 1 K = 1,3
          A(K+6) = XDOT(K,I)
          VI2 = VI2 + XDOT(K,I)**2
    1 CONTINUE
*
*       Include velocity kick in case of new WD, NS, BH or massless SN.
      IF (KW.GE.10.AND.KW.LE.15) THEN
          IKICK = .TRUE.
*       Ensure WD kicks are treated consistently (cf. define.f and #25).
          IF ((KW.EQ.10.OR.KW.EQ.11).AND.KZ(25).NE.1) IKICK = .FALSE.
          IF (KW.EQ.12.AND.KZ(25).NE.2) IKICK = .FALSE.
*       Distinguish between single star (first time only) and binary.
          IF (I.LE.N.AND.KW.NE.KSTAR(I).AND.IKICK) THEN
              CALL KICK(I,1,KW,DM)
          ELSE IF (I.GT.N.AND.IKICK) THEN
              IPAIR = I - N
              CALL KICK(IPAIR,0,KW,DM)
          END IF
      ELSE
          IKICK = .FALSE.
      END IF
*
*       Define consistent c.m. variables for KS mass loss (exclude Roche).
      IF (I.GT.N.AND.KSTAR(I).LE.10) THEN
          I2 = 2*(I - N)
          I1 = I2 - 1
          VF2 = 0.0
          DV2 = 0.0
          BODYI = BODY(I)
          IF (BODY(I).EQ.0.0D0) BODYI = BODY(I1) + BODY(I2)
          DO 8 K = 1,3
              X(K,I) = (BODY(I1)*X(K,I1) + BODY(I2)*X(K,I2))/BODYI
              XDOT(K,I) = (BODY(I1)*XDOT(K,I1) + BODY(I2)*XDOT(K,I2))/
     &                                                           BODYI
              X0(K,I) = X(K,I)
              X0DOT(K,I) = XDOT(K,I)
              VF2 = VF2 + XDOT(K,I)**2
              DV2 = DV2 + (XDOT(K,I) - A(K+6))**2
    8     CONTINUE
          VFAC = SQRT(VF2/VI2)
      END IF
*
*       Correct potential energy, all forces & first derivatives.
      POTJ = 0.0D0
      DO 40 J = IFIRST,NTOT
          IF (J.EQ.I) GO TO 40
          RIJ2 = 0.0D0
          RIJDOT = 0.0D0
          RDVDOT = 0.0D0
*
          DO 10 K = 1,3
              A(K) = X(K,I) - X(K,J)
              A(K+3) = A(K+6) - XDOT(K,J)
              RIJ2 = RIJ2 + A(K)**2
              RIJDOT = RIJDOT + A(K)*A(K+3)
              RDVDOT = RDVDOT + A(K)*(XDOT(K,I) - A(K+6))
   10     CONTINUE
*
          RIJ = SQRT(RIJ2)
          POTJ = POTJ + BODY(J)/RIJ
          A3 = 1.0/(RIJ2*RIJ)
          A4 = BODY(I)*A3
          A5 = DM*A3
          A6 = 3.0*RIJDOT/RIJ2
          A7 = 3.0*RDVDOT/RIJ2
*
          DO 15 K = 1,3
              A(K+3) = (A(K+3) - A6*A(K))*A5
              IF (IKICK) THEN
*       Include FDOT corrections due to increased velocity.
                  A(K+3) = A(K+3) + (XDOT(K,I) - A(K+6))*A4
                  A(K+3) = A(K+3) - A7*A(K)*A4
              END IF
   15     CONTINUE
*
*       Use neighbour list of #J to distinguish irregular & regular terms.
          NNB1 = LIST(1,J) + 1
          DO 25 L = 2,NNB1
              IF (LIST(L,J).EQ.I) THEN
                  DO 20 K = 1,3
                      F(K,J) = F(K,J) - 0.5*A(K)*A5
                      FI(K,J) = FI(K,J) - A(K)*A5
                      FDOT(K,J) = FDOT(K,J) - ONE6*A(K+3)
                      D1(K,J) = D1(K,J) - A(K+3)
                      FIDOT(K,J) = FIDOT(K,J) - A(K+3)
   20             CONTINUE
                  GO TO 40
              END IF
   25     CONTINUE      
          DO 30 K = 1,3
              F(K,J) = F(K,J) - 0.5*A(K)*A5
              FR(K,J) = FR(K,J) - A(K)*A5
              FDOT(K,J) = FDOT(K,J) - ONE6*A(K+3)
              D1R(K,J) = D1R(K,J) - A(K+3)
              FRDOT(K,J) = FRDOT(K,J) - A(K+3)
   30     CONTINUE
   40 CONTINUE
*
*       Update the potential and kinetic energy loss.
      EMDOT = EMDOT - DM*POTJ + 0.5*DM*VI2
*
*       Modify energy loss further for c.m. body (exclude Roche cases).
      IF (I.GT.N.AND.KSTAR(I).LE.10) THEN
          ECDOT = ECDOT - 0.5*BODY(I)*VI2*(VFAC**2 - 1.0)
      END IF
*
*       See whether linearized tidal terms should be included.
      IF (KZ(14).GT.0.AND.KZ(14).LT.3) THEN
          EMDOT = EMDOT - 0.5*DM*(TIDAL(1)*X(1,I)**2 +
     &                            TIDAL(3)*X(3,I)**2)
      END IF
*
*       Check optional Plummer potential.
      IF (KZ(14).EQ.4.OR.KZ(14).EQ.3) THEN
          RI2 = AP2
          DO 50 K = 1,3
              RI2 = RI2 + X(K,I)**2
  50      CONTINUE
          EMDOT = EMDOT - DM*MP/SQRT(RI2)
      END IF
*
*       Accumulate energy loss for conservation check (not used).
      E(12) = EMDOT
*
      RETURN
*
      END
back to top