https://github.com/florentrenaud/nbody6tt
Tip revision: 8aea50c213fd132d500c415511ae1e27eeabab80 authored by florent on 14 February 2015, 16:38:53 UTC
corrected mb typo in ttgalaxy
corrected mb typo in ttgalaxy
Tip revision: 8aea50c
tcirc.f
SUBROUTINE TCIRC(RP,ES0,I1,I2,ICIRC,TC)
*
*
* Circularization time.
* ---------------------
*
* Theory of Piet Hut, A & A 99, 126 (1981).
* Developed by Rosemary Mardling (31/1/97).
* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
INCLUDE 'common6.h'
REAL*8 WW(3),QQ(3),W(2),Q(2),AT0(2),M21,WG(2),QG(2),
& WSCALE(2),QSCALE(2),A(2),B(2),C(6)
DATA WW /2.119,3.113,8.175/
DATA QQ /0.4909,0.4219,0.2372/
DATA A /6.306505,-7.297806/
DATA B /32.17211,13.01598/
DATA C /5.101417,24.71539,-9.627739,1.733964,
& -2.314374,-4.127795/
DATA ECCM,ECCM1 /0.002,0.002001/
SAVE ITIME
DATA ITIME /0/
*
*
* Set large circularization time for merged binary or BHs.
IF (RADIUS(I1).EQ.0.0D0.OR.RADIUS(I2).EQ.0.0D0.OR.
& (KSTAR(I1).EQ.14.AND.KSTAR(I2).EQ.14)) THEN
TC = 1.0D+10
GO TO 30
END IF
*
* Specify index J1 as biggest radius to be used with AT0(1).
IF (RADIUS(I1).GE.RADIUS(I2)) THEN
J1 = I1
J2 = I2
ELSE
J1 = I2
J2 = I1
END IF
*
* Include an appropriate numerical factor (Chris Tout 07/08).
AR = RP/((1.0 - ES0)*RADIUS(J1))
FF = 1.0
IF (KSTAR(J1).EQ.1.AND.BODY(J1)*SMU.GT.1.25) THEN
FF = 100.0
ELSE IF (KSTAR(J1).EQ.4.OR.KSTAR(J1).EQ.7) THEN
FF = 100.0
ELSE IF (KSTAR(J1).GE.10) THEN
FF = 1.0D+17/AR**2
END IF
Q12 = BODY(J1)/BODY(J2)
TAUC = FF*2.0*Q12**2/(1.0 + Q12)*AR**8
* Exit on TC > 1.0D+08 yr (Chris Tout, Cambody Book 2008).
IF (TAUC.GT.1.0D+08) THEN
TC = 1.0D+10
GO TO 30
END IF
*
* Define oscillation period (dimensionless time) and damping constants.
XN = 0.0
DO 5 K = 1,2
IF (K.EQ.1) THEN
IK = J1
ELSE
IK = J2
END IF
* Specify polytropic index for each star (n = 3, 2 or 3/2).
IF (KSTAR(IK).EQ.3.OR.KSTAR(IK).EQ.5.OR.
& KSTAR(IK).EQ.6.OR.KSTAR(IK).EQ.9) THEN
IPAIR = KVEC(I1)
CALL GIANT(IPAIR,IK,WG,QG,WSCALE,QSCALE,XN,QL)
W(K) = WG(1)
Q(K) = QG(1)
ELSE
QL = 1.0D+04
IP = 3
IF (KSTAR(IK).GE.3) IP = 2
IF (KSTAR(IK).EQ.4.OR.KSTAR(IK).EQ.7) IP = 3
IF (KSTAR(IK).EQ.8) IP = 3
IF (KSTAR(IK).EQ.0) IP = 1
W(K) = WW(IP)
Q(K) = QQ(IP)
END IF
TL = TWOPI*RADIUS(IK)*SQRT(RADIUS(IK)/BODY(IK)/W(K))
AT0(K) = 1.0/(QL*TL)
5 CONTINUE
*
* Form mass, radius & pericentre ratio.
IF (RADIUS(I1).GE.RADIUS(I2)) THEN
M21 = BODY(I2)/BODY(I1)
R21 = RADIUS(I2)/RADIUS(I1)
RP1 = RP/RADIUS(I1)
ELSE
M21 = BODY(I1)/BODY(I2)
R21 = RADIUS(I1)/RADIUS(I2)
RP1 = RP/RADIUS(I2)
END IF
*
* Evaluate damping coefficient.
RR = RP1*(1.0 + ES0)
CONST = 2.0*(AT0(1)*(Q(1)/W(1))**2*(1.0 + M21)*M21 +
& AT0(2)*(Q(2)/W(2))**2*((1.0 + M21)/M21**2)*R21**8)/
& RR**8
*
* Adopt WD scaling for any NS/BH to avoid numerical problem.
IF (KSTAR(I1).GE.13.OR.KSTAR(I2).GE.13) THEN
CONST = 1.0D-04*CONST
END IF
*
* Form rational function approximation to Hut solution.
FF = (( A(2)*ES0 + A(1))*ES0 + 1.0 )/
& (( B(2)*ES0 + B(1))*ES0 + 1.0 )
* FF = MIN(FF,0.999D0)
*
* See whether we only want the modified eccentricity (routine BINPOP).
IF (ICIRC.LE.0) GO TO 10
*
TIME0 = 0.0
* Obtain the new eccentricity.
Z = (TIME - TIME0)*CONST + FF
ECC = (-1.0 + C(1)*Z - SQRT(C(2)*Z**2 + C(3)*Z + C(4)))
& /(C(5) + C(6)*Z)
*
ECC = MAX(ECC,ECCM)
RP = RP*(1.0 + ES0)/(1.0 + ECC)
ES0 = ECC
GO TO 30
*
* Evaluate circularization time (in units of 10**6 yrs).
10 TC = TSTAR*(1.0 - FF)/CONST
*
* Activate tidal indicator if TC < 2x10**9 yrs or hyperbolic orbit.
IF (TC.LT.2000.0.OR.ES0.GT.1.0) THEN
IP = KVEC(I1)
ITIME = ITIME + 1
* Note possibility of counter exceeding the limit.
IF (ITIME.GT.2000000000) ITIME = 0
IF (ICIRC.EQ.0.AND.KZ(27).EQ.2.AND.ITIME.LT.100) THEN
SEMI = -0.5*BODY(N+IP)/H(IP)
WRITE (6,20) I1, NCHAOS, ES0, RP1, M21, TC, SEMI, XN
20 FORMAT (' TCIRC: I1 NCH E RP M21 TC A n ',
& 2I5,F8.4,F8.1,F6.2,1P,2E10.2,0P,F5.1)
END IF
ICIRC = 1
* Define Roche search indicator for circularized orbit (ECCM1 > 0.002).
IF (ES0.LE.ECCM1.AND.KSTAR(N+IP).EQ.0) THEN
KSTAR(N+IP) = 10
END IF
ELSE
* Note ICIRC = -1 for some calls.
ICIRC = 0
END IF
*
30 RETURN
*
END