Revision dd8daa98413988b11c594e5a3bff8f5659d7ecf7 authored by Florent Renaud on 26 June 2014, 09:41:19 UTC, committed by Florent Renaud on 26 June 2014, 09:41:19 UTC
Originally released on 18 March 2013 Based on Nbody6 downloaded on 29 January 2013
1 parent edab1b7
bhstat.f
SUBROUTINE BHSTAT
*
*
* Black hole statistics.
* ----------------------
*
INCLUDE 'common6.h'
PARAMETER (NMX=10)
COMMON/CHREG/ TIMEC,TMAX,RMAXC,CX(10),NAMEC(NMX),NSTEP1,KZ27,KZ30
COMMON/CHAINC/ XC(3,NCMAX),UC(3,NCMAX),BODYC(NCMAX),ICH,
& LISTC(LMAX)
COMMON/POSTN2/ SEMIGR,ECCGR,DEGR,ISPIN
INTEGER NAMEB(10),IS(10)
REAL*8 CM(3),EI(10)
*
*
* Identify heavy single or binary in case of no chain.
IF (NCH.EQ.0) THEN
ICH = 0
DO 3 I = 1,NTOT
IF (BODY(I).GT.0.0001) ICH = I
3 CONTINUE
IF (ICH.EQ.0) GO TO 100
END IF
*
* Use c.m. as reference point unless case of strong recoil.
RICH = SQRT((X(1,ICH)-RDENS(1))**2 + (X(2,ICH)-RDENS(2))**2 +
& (X(3,ICH)-RDENS(3))**2)
IF (RICH.LT.0.2.AND.NCH.GT.0) THEN
DO 1 K = 1,3
CM(K) = X(K,ICH)
1 CONTINUE
ELSE
DO 2 K = 1,3
CM(K) = 0.0
2 CONTINUE
END IF
*
* Define names of all BHs.
Do 5 L = 1,NBH0
NAMEB(L) = L
5 CONTINUE
*
* Determine global location but skip ghosts.
LX = 0
DO 10 I = IFIRST,N
IF (BODY(I).EQ.0.0D0) GO TO 10
DO 8 L = 1,NBH0
IF (NAME(I).EQ.NAMEB(L)) THEN
LX = LX + 1
IS(LX) = I
IF (LX.EQ.NBH0) GO TO 12
GO TO 10
END IF
8 CONTINUE
10 CONTINUE
*
* Evaluate binding energy of at most LX single BHs and chain c.m.
12 IF (LX.EQ.0) GO TO 100
LX = LX + 1
* Ignore extra case if NCH = 0.
IS(LX) = ICH
DO 20 L = 1,LX
I = IS(L)
POTJ = 0.0D0
DO 16 J = IFIRST,NTOT
IF (J.EQ.I) GO TO 16
RIJ2 = 0.0D0
DO 15 K = 1,3
RIJ2 = RIJ2 + (X(K,I) - X(K,J))**2
15 CONTINUE
POTJ = POTJ + BODY(J)/SQRT(RIJ2)
16 CONTINUE
VI2 = 0.0
DO 18 K = 1,3
VI2 = VI2 + XDOT(K,I)**2
18 CONTINUE
* Note chain c.m. added for L = LX.
EI(L) = 0.5*VI2 - POTJ
IF (KZ(14).EQ.1) THEN
CALL XTRNLV(I,I)
EI(L) = EI(L) + HT/BODY(I)
END IF
20 CONTINUE
*
* Save central distance and binding energy for single BHs.
DO 30 L = 1,LX-1
I = IS(L)
RI = SQRT((X(1,I) - CM(1))**2 + (X(2,I) - CM(2))**2 +
& (X(3,I) - CM(3))**2)
RD = (X(1,I) - CM(1))*XDOT(1,I) + (X(2,I) - CM(2))*XDOT(2,I)
& + (X(3,I) - CM(3))*XDOT(3,I)
RD = RD/RI
RIJ2 = 0.0
VIJ2 = 0.0
RDOT = 0.0
DO 22 K = 1,3
RIJ2 = RIJ2 + (X(K,I) - X(K,ICH))**2
VIJ2 = VIJ2 + (XDOT(K,I) - XDOT(K,ICH))**2
RDOT = RDOT + (X(K,I) - X(K,ICH))*(XDOT(K,I)-XDOT(K,ICH))
22 CONTINUE
RIJ = SQRT(RIJ2)
SEMI = 2.0/RIJ - VIJ2/(BODY(I) + BODY(ICH))
SEMI = 1.0/SEMI
ECC2 = (1.0-RIJ/SEMI)**2 + RDOT**2/(SEMI*(BODY(I)+BODY(ICH)))
ECC = SQRT(ECC2)
IF (RI.GT.100.0*RMIN) ECC = 0.0
PMIN = SEMI*(1.0 - ECC)/SEMIGR
IF (RI.GT.100.0*RMIN) PMIN = 0.0
IF (NAME(I).EQ.1) THEN
WRITE (91,25) TTOT, ECC, PMIN, RI, EI(L)
25 FORMAT (' ',F8.2,F7.3,F6.1,1P,2E10.2)
CALL FLUSH(91)
ELSE IF (NAME(I).EQ.2) THEN
WRITE (92,25) TTOT, ECC, PMIN, RI, EI(L)
CALL FLUSH(92)
ELSE IF (NAME(I).EQ.3) THEN
WRITE (93,25) TTOT, ECC, PMIN, RI, EI(L)
CALL FLUSH(93)
ELSE IF (NAME(I).EQ.4) THEN
WRITE (94,25) TTOT, ECC, PMIN, RI, EI(L)
CALL FLUSH(94)
END IF
30 CONTINUE
*
* Assign small values inside the chain for plotting purposes (log!).
ECC = 0.0
PMIN = 0.0
DO 40 L = 1,NCH
RI = 1.0D-04
IF (NAMEC(L).EQ.1) THEN
WRITE (91,25) TTOT, ECC, PMIN, RI, EI(LX)
CALL FLUSH(91)
ELSE IF (NAMEC(L).EQ.2) THEN
WRITE (92,25) TTOT, ECC, PMIN, RI, EI(LX)
CALL FLUSH(92)
ELSE IF (NAMEC(L).EQ.3) THEN
WRITE (93,25) TTOT, ECC, PMIN, RI, EI(LX)
CALL FLUSH(93)
ELSE IF (NAMEC(L).EQ.4) THEN
WRITE (94,25) TTOT, ECC, PMIN, RI, EI(LX)
CALL FLUSH(94)
END IF
40 CONTINUE
*
* Include data for any chain c.m.
IF (NCH.GT.0) THEN
VI2 = XDOT(1,ICH)**2 + XDOT(2,ICH)**2 + XDOT(3,ICH)**2
WRITE (95,50) TTOT, ECC, PMIN, RICH, EI(LX), VI2
50 FORMAT (' ',F8.2,F7.3,F6.1,1P,2E10.2,0P,F7.3)
CALL FLUSH(95)
END IF
*
100 RETURN
*
END
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...