https://github.com/nbodyx/Nbody6
Tip revision: dff3b5c68673d4d4c4c9f54c8ba110b8416098f2 authored by nitadori on 31 May 2020, 13:04:00 UTC
Avoided touching uninitialized SEMIX
Avoided touching uninitialized SEMIX
Tip revision: dff3b5c
ficorr.f
SUBROUTINE FICORR(I,DM)
*
*
* Local force corrections due to mass loss.
* -----------------------------------------
*
INCLUDE 'common6.h'
REAL*8 A(6)
*
*
* Correct force & first derivative for neighbours of body #I.
NNB1 = LIST(1,I) + 1
DO 20 L = 2,NNB1
J = LIST(L,I)
NNB2 = LIST(1,J) + 1
*
* Ensure that body #I is a neighbour of body #J.
DO 1 K = 2,NNB2
IF (LIST(K,J).EQ.I) GO TO 5
1 CONTINUE
GO TO 20
*
5 RIJ2 = 0.0D0
A7 = 0.0D0
*
DO 10 K = 1,3
A(K) = X(K,I) - X(K,J)
A(K+3) = XDOT(K,I) - XDOT(K,J)
RIJ2 = RIJ2 + A(K)**2
A7 = A7 + A(K)*A(K+3)
10 CONTINUE
*
A7 = 3.0*A7/RIJ2
A8 = DM/(RIJ2*SQRT(RIJ2))
*
DO 15 K = 1,3
A(K+3) = (A(K+3) - A7*A(K))*A8
F(K,J) = F(K,J) - 0.5*A(K)*A8
FI(K,J) = FI(K,J) - A(K)*A8
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)
15 CONTINUE
20 CONTINUE
*
* Obtain the potential energy due to all particles.
POTJ = 0.0D0
DO 30 J = IFIRST,NTOT
IF (J.EQ.I) GO TO 30
RIJ2 = 0.0D0
DO 25 K = 1,3
RIJ2 = RIJ2 + (X(K,I) - X(K,J))**2
25 CONTINUE
POTJ = POTJ + BODY(J)/SQRT(RIJ2)
30 CONTINUE
*
* Update the potential energy loss.
EMDOT = EMDOT - DM*POTJ
*
* Include kinetic energy correction.
VI2 = XDOT(1,I)**2 + XDOT(2,I)**2 + XDOT(3,I)**2
EMDOT = EMDOT + 0.5*DM*VI2
*
* See whether linearized tidal terms should be included.
IF (KZ(14).GT.0.AND.KZ(14).LT.3) THEN
ECDOT = ECDOT - 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
ECDOT = ECDOT - DM*MP/SQRT(RI2)
END IF
*
RETURN
*
END