https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
hcorr.f
SUBROUTINE HCORR(I,DM,RNEW)
*
*
* Mass loss correction of KS orbit.
* ---------------------------------
*
INCLUDE 'common6.h'
COMMON/BINARY/ CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
& HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
& NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAG(MMAX)
*
*
* Copy pair index and determine secondary component (note I1 = I here).
IPAIR = KSPAIR
I1 = I
I2 = 2*IPAIR - 1
IF (I2.EQ.I1) THEN
I2 = I2 + 1
END IF
*
* Define c.m. index and old semi-major axis.
J = N + IPAIR
SEMI0 = -0.5D0*BODY(J)/H(IPAIR)
*
* Distinguish different cases (KS, triple or outer ghost of quadruple).
IF (BODY(I).GT.0.0D0) THEN
ZMU0 = BODY(I1)*BODY(I2)/BODY(J)
ELSE
IM = 0
JJ = I
IF (I.LT.IFIRST) JJ = N + KVEC(I)
* Determine merger index for inner and outer binary masses.
DO 1 K = 1,NMERGE
IF (NAMEG(K).EQ.NAME(JJ)) IM = K
1 CONTINUE
IF (IM.EQ.0) THEN
WRITE (6,5) NAME(I), NAME(JJ), (NAMEG(K),K=1,NMERGE)
5 FORMAT (' DANGER HCORR! NMI NMJJ NMG ',12I6)
STOP
END IF
* Form old reduced mass from inner and outer binary.
ZM1 = CM(1,IM) + CM(2,IM)
ZM2 = CM(3,IM) + CM(4,IM)
ZMU0 = ZM1*ZM2/BODY(J)
END IF
*
* Obtain energy change from MK*A = const and H = -M/(2*A) (DCH 8/96).
DH = DM/SEMI0*(1.0 - 0.5*DM/BODY(J))
*
* Reduce mass of c.m. body and subtract energy correction terms.
BODY(J) = BODY(J) - DM
IF (BODY(I).GT.0.0D0) THEN
ZMU1 = (BODY(I1) - DM)*BODY(I2)/BODY(J)
ELSE
ZMU1 = ZM1*(ZM2 - DM)/BODY(J)
END IF
EMDOT = EMDOT - ZMU1*DH - (ZMU1 - ZMU0)*H(IPAIR)
* Retain binding energy gained in past interactions.
EGRAV = EGRAV - ZMU1*DH - (ZMU1 - ZMU0)*H(IPAIR)
*
* Include diagnostic warning for large relative energy change.
IF (DH.GT.0.2*ABS(H(IPAIR))) THEN
WRITE (6,10) NAME(I), DH, H(IPAIR), R(IPAIR)/SEMI0, DM*SMU
10 FORMAT (' WARNING! LARGE CORRECTION NM DH H R/A DMS ',
& I6,2F8.2,2F7.2)
END IF
*
* Update the binding energy due to mass loss DM and set new radius.
H(IPAIR) = H(IPAIR) + DH
RADIUS(I) = RNEW
*
* Modify KS variables at constant eccentricity.
CALL EXPAND(IPAIR,SEMI0)
*
RETURN
*
END