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
adjust.f
SUBROUTINE ADJUST
*
*
* Parameter adjustment and energy check.
* --------------------------------------
*
INCLUDE 'common6.h'
COMMON/ECHAIN/ ECH
SAVE DTOFF
DATA DTOFF /100.0D0/
*
*
* Predict X & XDOT for all particles (except unperturbed pairs).
CALL XVPRED(IFIRST,NTOT)
*
* Obtain the total energy at current time (resolve all KS pairs).
CALL ENERGY
*
* Initialize c.m. terms.
DO 10 K = 1,3
CMR(K) = 0.0D0
CMRDOT(K) = 0.0D0
10 CONTINUE
*
* Obtain c.m. & angular momentum integrals and Z-moment of inertia.
AZ = 0.0D0
ZM = 0.0D0
ZMASS = 0.0D0
ISUN = 1
DO 20 I = 1,N
IF (NAME(I).EQ.1) ISUN = I
ZMASS = ZMASS + BODY(I)
DO 15 K = 1,3
CMR(K) = CMR(K) + BODY(I)*X(K,I)
CMRDOT(K) = CMRDOT(K) + BODY(I)*XDOT(K,I)
15 CONTINUE
* RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
* & (X(3,I) - RDENS(3))**2
AZ = AZ + BODY(I)*(X(1,I)*XDOT(2,I) - X(2,I)*XDOT(1,I))
ZM = ZM + BODY(I)*(X(1,I)**2 + X(2,I)**2)
20 CONTINUE
*
* Form c.m. coordinates & velocities (vectors & scalars).
DO 25 K = 1,3
CMR(K) = CMR(K)/ZMASS
CMRDOT(K) = CMRDOT(K)/ZMASS
25 CONTINUE
*
CMR(4) = SQRT(CMR(1)**2 + CMR(2)**2 + CMR(3)**2)
CMRDOT(4) = SQRT(CMRDOT(1)**2 + CMRDOT(2)**2 + CMRDOT(3)**2)
*
* Subtract the kinetic energy of c.m. due to possible cloud effects.
IF (KZ(13).GT.0) ZKIN = ZKIN - 0.5*ZMASS*CMRDOT(4)**2
*
* Include non-zero virial energy for Plummer potential and/or 3D case.
VIR = POT - VIR
* Note angular momentum term is contained in virial energy (#14=1/2).
Q = ZKIN/VIR
E(3) = ZKIN - POT + EPL
* Modify single particle energy by tidal energy (except pure 3D).
IF (KZ(14).NE.3) THEN
E(3) = E(3) + ETIDE
END IF
*
* Modify angular momentum integral using Chandrasekhar eq. (5.530).
IF (KZ(14).EQ.1.OR.KZ(14).EQ.2) THEN
AZ = AZ + 0.5*TIDAL(4)*ZM
END IF
*
* Define crossing time using single particle energy (cf. option 14).
TCR = ZMASS**2.5/(2.0*ABS(E(3)))**1.5
* Note: see below for better definition in Plummer or 3D orbit case.
IF (Q.GT.1.0.AND.KZ(14).LT.3) THEN
TCR = TCR*SQRT(2.0*Q)
END IF
* Form provisional total energy.
ETOT = ZKIN - POT + ETIDE + EPL
*
* Include KS pairs, triple, quad, mergers, collisions & chain.
ETOT = ETOT + EBIN + ESUB + EMERGE + ECOLL + EMDOT + ECDOT
* Add any contributions from chain regularization.
IF (NCH.GT.0) THEN
ETOT = ETOT + ECH
END IF
*
* Update energies and form the relative error (divide by ZKIN or E(3)).
IF (TIME.LE.0.0D0) THEN
DE = 0.0D0
BE(1) = ETOT
BE(3) = ETOT
DELTA1 = 0.0D0
ELSE
BE(2) = BE(3)
BE(3) = ETOT
DE = BE(3) - BE(2)
DELTA1 = DE
DETOT = DETOT + DE
DE = DE/MAX(ZKIN,ABS(E(3)))
* Save sum of relative energy error for main output and accumulate DE.
ERROR = ERROR + DE
ERRTOT = ERRTOT + DE
END IF
*
* Set provisional half-mass radius.
RSCALE = 0.5*ZMASS**2/POT
*
* Determine average neighbour number and smallest neighbour sphere.
NNB = 0
RS0 = RSCALE
DO 30 I = IFIRST,NTOT
NNB = NNB + LIST(1,I)
IF (LIST(1,I).GT.0) RS0 = MIN(RS0,RS(I))
30 CONTINUE
NNB = NNB/(N - NPAIRS)
*
* Use current value if minimum neighbour sphere not implemented.
IF (RSMIN.EQ.0.0D0) RSMIN = RS0
*
* Find density centre & core radius (Casertano & Hut, Ap.J. 298, 80).
IF (N-NPAIRS.GE.20.AND.KZ(29).EQ.0.AND.KZ(5).NE.3) THEN
CALL CORE
ELSE
NC = N
ZMC = ZMASS
RHOD = 1.0
RHOM = 1.0
RC = RSCALE
RC22 = RC**2
RC2IN = 1.0/RC22
END IF
*
* Take the Sun as reference for plotting planetesimal disk members.
IF (KZ(5).EQ.3) THEN
DO 35 K = 1,3
RDENS(K) = X(K,ISUN)
35 CONTINUE
*
* Determine the eccentricity dispersion and total energy of disk.
DISP2 = 0.0
EDISK = 0.0
DO 40 I = 1,N
IF (NAME(I).EQ.1.OR.NAME(I).EQ.NZERO) GO TO 40
RI2 = (X(1,I) - X(1,ISUN))**2 + (X(2,I) - X(2,ISUN))**2
VI2 = (XDOT(1,I) - XDOT(1,ISUN))**2 +
& (XDOT(2,I) - XDOT(2,ISUN))**2
RRDOT = (X(1,I) - X(1,ISUN))*(XDOT(1,I) - XDOT(1,ISUN)) +
& (X(2,I) - X(2,ISUN))*(XDOT(2,I) - XDOT(2,ISUN))
RI = SQRT(RI2)
SEMI = 2.0/RI - VI2/(BODY(ISUN) + BODY(I))
SEMI = 1.0/SEMI
ECC2 = (1.0 - RI/SEMI)**2 +
& RRDOT**2/(SEMI*(BODY(I) + BODY(ISUN)))
DISP2 = DISP2 + ECC2
EDISK = EDISK - 0.5*BODY(I)/SEMI
40 CONTINUE
DISP = SQRT(DISP2/FLOAT(N-2))
WRITE (35,42) TTOT, DISP, EDISK
42 FORMAT (' ',F8.1,1P,E10.2,E12.4)
END IF
*
* Check optional sorting of Lagrangian radii & half-mass radius.
IF (KZ(7).GT.0) THEN
CALL LAGR(RDENS)
END IF
*
* Scale average & maximum core density by the mean value.
RHOD = 4.0*TWOPI*RHOD*RSCALE**3/(3.0*ZMASS)
RHOM = 4.0*TWOPI*RHOM*RSCALE**3/(3.0*ZMASS)
*
* Adopt density contrasts of unity for hot system.
IF (KZ(29).GT.0.AND.ZKIN.GT.POT) THEN
RHOD = 1.0
RHOM = 1.0
END IF
*
* Check optional determination of regularization parameters.
IF (KZ(16).GT.0) THEN
RMIN0 = RMIN
*
* Form close encounter distance from scale factor & density contrast.
RMIN = 4.0*RSCALE/(FLOAT(N)*RHOD**0.3333)
* Include alternative expression based on core radius (experimental).
IF (KZ(16).GT.1.AND.NC.LT.0.01*N) THEN
RMIN = 0.01*RC/FLOAT(NC)**0.3333
END IF
* Use harmonic mean to reduce fluctuations (avoid initial value).
IF (TIME.GT.0.0D0) RMIN = SQRT(RMIN0*RMIN)
* Impose maximum value for sufficient perturbers.
RMIN = MIN(RMIN,RSMIN*GMIN**0.3333)
* Define scaled DTMIN by RMIN & <M> and include ETAI for consistency.
DTMIN = 0.04*SQRT(ETAI/0.02D0)*SQRT(RMIN**3/BODYM)
* Specify binding energy per unit mass of hard binary (impose Q = 0.5).
ECLOSE = 4.0*MAX(ZKIN,ABS(ZKIN - POT))/ZMASS
* Adopt central velocity as upper limit (avoids large kick velocities).
IF (2.0*ZKIN/ZMASS.GT.VC**2.AND.VC.GT.0.0) ECLOSE = 2.0*VC**2
IF (Q.GT.0.5) THEN
ECLOSE = 0.5*ECLOSE/Q
END IF
ECLOSE = MIN(ECLOSE,1.0D0)
KSMAG = 0
END IF
*
* Check optional modification of DTMIN, ECLOSE & TCR for hot system.
IF (KZ(29).GT.0.AND.Q.GT.1.0) THEN
DTMIN = 0.04*SQRT(ETAI/0.02D0)*SQRT(RMIN**3/BODYM)
SIGMA2 = 2.0*ZKIN/ZMASS
VP2 = 4.0*BODYM/RMIN
DTMIN = DTMIN*SQRT((VP2 + SIGMA2/Q)/(VP2 + 2.0D0*SIGMA2))
ECLOSE = SIGMA2
TCR = 2.0*RSCALE/SQRT(SIGMA2)
END IF
*
* Set useful scalars for the integrator.
SMIN = 2.0*DTMIN
RMIN2 = RMIN**2
RMIN22 = 4.0*RMIN2
EBH = -0.5*BODYM*ECLOSE
IF (TIME.LE.0.0D0) THEN
STEPJ = 0.01*(60000.0/FLOAT(N))**0.3333
IF (DMOD(SMAX,DTK(10)).NE.0.0D0) THEN
WRITE (6,43) SMAX, SMAX/DTK(10)
43 FORMAT (' FATAL ERROR! SMAX SMAX/DTK(10) ',1P,2E10.2)
STOP
END IF
END IF
* Adopt 2*RSMIN for neighbour sphere volume factor in routine REGINT.
RSFAC = MAX(25.0/TCR,3.0D0*VC/(2.0D0*RSMIN))
*
* Update density contrast factor for neighbour sphere modification.
IF (TIME.LE.0.0D0.OR.KZ(40).EQ.0) THEN
ALPHA = FLOAT(NNBMAX)*SQRT(0.08D0*RSCALE**3/FLOAT(N-NPAIRS))
END IF
* Include optional stabilization to increase neighbour number.
IF (KZ(40).EQ.1.AND.FLOAT(NNB).LT.0.5*NNBMAX) THEN
FAC = 1.0 + (0.5*NNBMAX - NNB)/(0.5*FLOAT(NNBMAX))
ALPHA = FAC*ALPHA
END IF
*
* Define tidal radius for isolated system (2*RTIDE used in ESCAPE).
IF (KZ(14).EQ.0) RTIDE = 10.0*RSCALE
*** FlorentR - set the tidal radius
IF (KZ(14).EQ.9) RTIDE = 10.0*RSCALE
*** FRenaud
* Redefine the crossing time for 3D cluster orbit or Plummer model.
IF ((KZ(14).EQ.3.OR.KZ(14).EQ.4).AND.ZKIN.GT.0.0) THEN
TCR = 2.0*RSCALE/SQRT(2.0*ZKIN/ZMASS)
END IF
*
* Print energy diagnostics & KS parameters.
ICR = TTOT/TCR
WRITE (6,45) TTOT, Q, DE, BE(3), RMIN, DTMIN, ICR, DELTA1, E(3),
& DETOT
45 FORMAT (/,' ADJUST: TIME =',F8.2,' Q =',F5.2,' DE =',1P,E9.1,
& ' E =',0P,F10.6,' RMIN =',1P,E8.1,' DTMIN =',E8.1,
& ' TC =',0P,I5,' DELTA =',1P,E9.1,' E(3) =',0P,F10.6,
& ' DETOT =',F10.6)
CALL FLUSH(6)
*
* Perform automatic error control (RETURN on restart with KZ(2) > 1).
CALL CHECK(DE)
IF (ABS(DE).GT.5.0*QE) GO TO 70
*
* Check for escaper removal.
IF (KZ(23).GT.0) THEN
CALL ESCAPE
END IF
*
* Check correction for c.m. displacements.
IF (KZ(31).GT.0) THEN
CALL CMCORR
END IF
*
* Include diagnostics for massive binary (bound or unbound initially).
IF (KZ(5).EQ.4) THEN
IP = 0
DO 50 IPAIR = 1,NPAIRS
IF (NAME(2*IPAIR-1).LE.2.OR.NAME(2*IPAIR).LE.2) THEN
IP = IPAIR
END IF
50 CONTINUE
IF (IP.GT.0) THEN
I1 = 2*IP - 1
I2 = I1 + 1
SEMI = -0.5*BODY(N+IP)/H(IP)
ECC2 = (1.0 - R(IP)/SEMI)**2 +
& TDOT2(IP)**2/(SEMI*BODY(N+IP))
EB = BODY(I1)*BODY(I2)*H(IP)/BODY(N+IP)
WRITE (35,52) TTOT, SEMI, EB, E(3), SQRT(ECC2),
& NAME(I1), NAME(I2)
52 FORMAT (' ',F8.1,1P,3E12.4,0P,F7.3,2I5)
* 52 FORMAT (' T A E EB ECL NAME ',F8.1,1P,3E12.4,0P,F7.3,2I5)
CALL FLUSH(35)
END IF
END IF
*
* Allow a gradual decrease of NNBMAX due to escaper removal.
IF (KZ(40).EQ.3) THEN
NNBMAX = NBZERO*SQRT(FLOAT(N)/FLOAT(NZERO))
ZNBMAX = 0.9*NNBMAX
ZNBMIN = MAX(0.2*NNBMAX,1.0)
END IF
*
* Include optional fine-tuning of neighbour number (#40 >= 2).
IF (KZ(40).GE.2) THEN
IF (NNB.LT.NNBMAX/5.0) THEN
ALPHA = 1.1*ALPHA
ELSE
ALPHA = 0.9*ALPHA
END IF
END IF
*
* See whether standard output is due.
IOUT = 0
IF (TIME.GE.TNEXT) THEN
CALL OUTPUT
IOUT = 1
END IF
*
* Include optional diagnostics for the hardest binary below ECLOSE.
IF (KZ(33).GE.2.AND.IOUT.GT.0) THEN
HP = 0.0
IP = 0
DO 60 IPAIR = 1,NPAIRS
* Skip outer (ghost) binary of quadruple system.
IF (H(IPAIR).LT.HP.AND.BODY(N+IPAIR).GT.0.0D0) THEN
HP = H(IPAIR)
IP = IPAIR
END IF
60 CONTINUE
IF (IP.GT.0.AND.HP.LT.-ECLOSE) THEN
I1 = 2*IP - 1
I2 = I1 + 1
SEMI = -0.5*BODY(N+IP)/H(IP)
PB = DAYS*SEMI*SQRT(SEMI/BODY(N+IP))
ECC2 = (1.0 - R(IP)/SEMI)**2 +
& TDOT2(IP)**2/(SEMI*BODY(N+IP))
EB = BODY(I1)*BODY(I2)*H(IP)/BODY(N+IP)
WRITE (39,62) TTOT, NAME(I1), NAME(I2), KSTAR(N+IP),
& LIST(1,I1), SQRT(ECC2), SEMI, PB, EB, E(3)
62 FORMAT (' BINARY: T NAME K* NP E A P EB E3 ',
& F8.1,2I6,2I4,F7.3,1P,2E10.2,0P,2F9.4)
CALL FLUSH(39)
END IF
END IF
*
* Include optional KS reg of binaries outside standard criterion.
IF (KZ(8).GT.0.AND.DMOD(TIME,DTK(10)).EQ.0.0D0) THEN
* Note DMOD condition needed for CALL KSREG and CALL STEPS.
DTCL = 30.0*DTMIN
RCL = 10.0*RMIN
CALL SWEEP(DTCL,RCL)
END IF
*
* Update time for next adjustment.
TADJ = TADJ + DTADJ
* Re-initialize marginal stability counter to avoid including old case.
NMTRY = 0
*
* Check optional truncation of time.
IF (KZ(35).GT.0.AND.TIME.GE.DTOFF) THEN
CALL OFFSET(DTOFF)
END IF
*
* Obtain elapsed CPU time and update total since last output/restart.
CALL CPUTIM(TCOMP)
CPUTOT = CPUTOT + TCOMP - CPU0
CPU0 = TCOMP
*
* Save COMMON after energy check (skip TRIPLE, QUAD, CHAIN).
TDUMP = TIME
IF (KZ(2).GE.1.AND.NSUB.EQ.0) CALL MYDUMP(1,2)
*
* Check termination criteria (TIME > TCRIT & N <= NCRIT).
IF (TTOT.GE.TCRIT.OR.N.LE.NCRIT.OR.TTOT+DTADJ.GT.TCRIT) THEN
* Terminate after optional COMMON save.
WRITE (6,65) TTOT, CPUTOT/60.0, ERRTOT, DETOT
65 FORMAT (//,9X,'END RUN',3X,'TIME =',F7.1,' CPUTOT =',F7.1,
& ' ERRTOT =',F10.6,' DETOT =',F10.6)
IF (KZ(1).GT.0.AND.NSUB.EQ.0) CALL MYDUMP(1,1)
STOP
END IF
*
70 RETURN
*
END