https://github.com/florentrenaud/nbody6tt
Raw File
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
chrect.f
      SUBROUTINE CHRECT(IPAIR,DMR)
*
*
*       Rectification of chaotic orbit.
*       -------------------------------
*
      INCLUDE 'common6.h'
      COMMON/MODES/  EB0(NTMAX),ZJ0(NTMAX),ECRIT(NTMAX),AR(NTMAX),
     &               BR(NTMAX),EOSC(4,NTMAX),EDEC(NTMAX),TOSC(NTMAX),
     &               RP(NTMAX),ES(NTMAX),CM(2,NTMAX),IOSC(NTMAX),
     &               NAMEC(NTMAX)
      REAL*8  WW(6),W(4),WG(2),QG(2),WSCALE(2),QSCALE(2)
      CHARACTER*8  WHICH1
      LOGICAL  SLEEP
      DATA  WW  /2.119,3.113,8.175,3.742,4.953,9.413/
*
*
*       Define c.m. & KS indices and search current names for chaos index.
      I = N + IPAIR
      I1 = 2*IPAIR - 1
      I2 = I1 + 1
      IC = 0
      DO 1 K = 1,NCHAOS
          IF (NAMEC(K).EQ.NAME(I)) IC = K
    1 CONTINUE
*
*       Include case of chain chaos without identified NAMEC.
      IF (IC.EQ.0.AND.NCHAOS.GT.0) THEN
          WRITE (6,2)  NCHAOS, IPAIR, KSTAR(I), NAME(I1), NAME(I2),
     &                 LIST(1,I1), STEP(I1), STEP(I)
    2     FORMAT (' WARNING!    CHRECT    NCH KS K* NAM NP DT1 DTI ',
     &                                    3I4,2I6,I4,1P,2E10.2)
*       See whether wrong component name (+ NZERO) saved as NAMEC in CHAOS2.
          NAM2 = 0
          DO 6 L = 1,2
              NAM1 = KSAVE(2*L) - NAME(I2)
              IF (KSAVE(2*L-1).LT.0.AND.NAM1.EQ.NAME(I1)) THEN
                  NAM2 = NZERO + NAME(I2)
              END IF
    6     CONTINUE
          IC = NCHAOS
          NAMC = NAMEC(IC)
*       Check identification for correct value (two K*=-2 are possible).
          DO 8 K = 1,NCHAOS
              IF (NAMEC(K).EQ.NAM2) IC = K
    8     CONTINUE
          NAMEC(IC) = NAME(I)
          IF (NAM2.EQ.NAMC.OR.IC.LT.NCHAOS) THEN
              WRITE (6,9)  IC, NCHAOS, NAM2, NAMC, NAME(I)
    9         FORMAT (' CHRECT RESTORE    IC NCH NM2 NMC NMI ',2I4,3I8)
          END IF
      ELSE IF (NCHAOS.EQ.0) THEN
          WRITE (6,3)  NCHAOS, IPAIR, KSTAR(I), NAME(I)
    3     FORMAT (' CHRECT RESTORE    NCH KS K* NAM ',3I4,I6)
*       Restore case of former merger with KSTARM < 0 to chaos table.
          NCHAOS = 1
          IC = 1
          NAMEC(NCHAOS) = NAME(I)
      END IF
*
*       Save variables for diagnostic output.
      TIME0 = TOSC(IC)
      ES0 = ES(IC)
      TC = 0.0
      IDIS = KSTAR(I)
*
*       Obtain current values of KS variables in case of spiral.
      IF (KSTAR(I).EQ.-2) THEN
*       Skip on call from RESET/MDOT with arbitrary phase (only updating).
          IF (ABS(TDOT2(IPAIR)).GT.1.0D-12.AND.DMR.GE.0.0) THEN
*       Rectify KS variables in order to obtain correct pericentre.
              CALL KSRECT(IPAIR)
*       Reduce eccentric anomaly by pi for inward motion.
              IF (TDOT2(IPAIR).LT.0.0D0) THEN
                  CALL KSAPO(IPAIR)
              END IF
*       Transform from outward motion to exact pericentre.
              CALL KSPERI(IPAIR)
          END IF
*
*       Form current two-body elements.
          SEMI = -0.5*BODY(I)/H(IPAIR)
          ECC2 = (1.0 - R(IPAIR)/SEMI)**2 +
     &                                    TDOT2(IPAIR)**2/(BODY(I)*SEMI)
          ECC = SQRT(ECC2)
*
*       Update periastron and eccentricity (ECC modulation or mass loss).
          QPERI = SEMI*(1.0D0 - ECC)
          RP(IC) = QPERI
          ES(IC) = ECC
*
*       Update orbital parameters after merger, mass loss or radius change.
          IF (DMR.GE.0.0D0) THEN
              KSTAR(I) = -KSTAR(I)
              CALL SPIRAL(IPAIR)
              IF (IPHASE.LT.0.OR.KSTAR(I).GT.0) GO TO 30
          END IF
*
*       Check circularization time for spiral after radius/orbit expansion.
          ICIRC = -1
          CALL TCIRC(QPERI,ECC,I1,I2,ICIRC,TC)
          TC = MAX(TC,0.01D0)
*       Restrict lookup time to TC/2 for exit from possible merger.
          TEV(I1) = MIN(TEV(I1),TIME + 0.5*TC/TSTAR)
*         TEV(I2) = TEV(I1)
      END IF
*
*       Form (new) semi-major axis, eccentricity & periastron distance.
      SEMI = -0.5*BODY(I)/H(IPAIR)
      ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
      ECC = SQRT(ECC2)
      QPERI = SEMI*(1.0D0 - ECC)
*
*       Set latest values of periastron & eccentricity and update time.
      RP(IC) = QPERI
      ES(IC) = ECC
      TOSC(IC) = TIME
*
*       Include check for SLEEP after recent WD/NS formation.
      SLEEP = .FALSE.
      KM = MAX(KSTAR(I1),KSTAR(I2))
      IF (KM.GE.10) THEN
*       Determine index of most recent degenerate object formation.
          IF (KSTAR(I2).LT.10) THEN
              J1 = I1
          ELSE IF (KSTAR(I1).LT.10) THEN
              J1 = I2
          ELSE
              J1 = I1
              IF (EPOCH(I2).GT.EPOCH(I1)) J1 = I2
          END IF
          IF (TIME - TEV0(J1).LT.2.0*STEPX.AND.TC.GT.3000.0) THEN
              SLEEP = .TRUE.
          END IF
      END IF
*
*       Check circularization time after chain regularization.
      IF (IPHASE.EQ.8) THEN
          ICIRC = -1
          CALL TCIRC(QPERI,ECC,I1,I2,ICIRC,TC)
          IF (TC.GT.3000.0) THEN
              WRITE (6,4)  NAME(I1), TC
    4         FORMAT (' CHAIN SLEEP:    NM TC ',I6,1P,E10.2)
              SLEEP = .TRUE.
          END IF
      END IF
*
*       Re-initialize all chaos parameters after expansion or check SLEEP.
      IF (DMR.GT.0.01.OR.KSTAR(I).EQ.-1.OR.SLEEP) THEN
          IF (DMR.GT.0.01.AND.TC.GT.100.0.AND.KM.LT.5) THEN
              WRITE (6,5)  NAME(I1), NAME(I2), KSTAR(I1), KSTAR(I2),
     &                     KSTAR(I), TPHYS, RADIUS(I1), RADIUS(I2),
     &                     QPERI, SEMI, ECC, ES0, BODY(I)*ZMBAR, TC
    5         FORMAT (' CHRECT:    NAM K* TP R* QP A E E0 M TC ',
     &                             2I6,3I4,F8.1,1P,4E10.2,0P,3F7.3,F7.1)
          END IF
*
*       Reset spiral indicator for degenerate component and long t_{circ}.
          IF (SLEEP) THEN
              KSTAR(I) = 0
              NSLP = NSLP + 1
              TK = SEMI*SQRT(SEMI/BODY(I))
              TB = YRS*TK
              XP = (TIME - TIME0)/(TWOPI*TK)
              QPS = SEMI*(1.0 - ECC)/MAX(RADIUS(I1),RADIUS(I2))
              WRITE (6,10)  TTOT, NAME(I1), NAME(I2), KSTAR(I1), 
     &                      KSTAR(I2), ECC, ES0, QPS, SEMI, TC, TB, XP
   10         FORMAT (' SLEEP SPIRAL    T NM K* E E0 QP/S A TC TB DTK ',
     &                                  F9.2,2I6,2I4,2F8.4,1P,5E9.1)
*       Update chaos variables at end of routine SPIRAL (argument < 0).
              II = -I
              CALL SPIRAL(II)
              GO TO 30
          END IF
*
*       Save eccentricity, binding energy & J0 and initialize EDEC & IOSC.
          ZMU = BODY(I1)*BODY(I2)/BODY(I)
          CJ = ZMU*SQRT(BODY(I))
          EB0(IC) = ZMU*H(IPAIR)
          ZJ0(IC) = CJ*SQRT(QPERI*(1.0 + ECC))
          EDEC(IC) = 0.0
          IOSC(IC) = 1
          XN = 0.0
          DO 12 K = 1,2
              IK = I1 + K - 1
              IF (KSTAR(IK).EQ.3.OR.KSTAR(IK).EQ.5.OR.
     &            KSTAR(IK).EQ.6.OR.KSTAR(IK).EQ.9) THEN
                  CALL GIANT(IPAIR,IK,WG,QG,WSCALE,QSCALE,XN,QL)
                  W(K) = WG(1)
              ELSE
                  IP = 3
                  IF (KSTAR(IK).EQ.0) IP = 1
                  W(K) = WW(IP)
              END IF
   12     CONTINUE
*
*       Set new chaos boundary parameters (ECRIT, AR & BR).
          CALL CHAOS0(QPERI,ECC,EB0(IC),ZJ0(IC),BODY(I1),BODY(I2),
     &             RADIUS(I1),RADIUS(I2),W,ECRIT(IC),AR(IC),BR(IC),IDIS)
*
          RCOLL = RADIUS(I1) + RADIUS(I2)
	  IF (IDIS.EQ.-1.AND.KSTAR(I).EQ.-1) THEN
              IOSC(IC) = 2
              WRITE (6,15)  TTOT, IPAIR, NAME(I1), NAME(I2), KSTAR(I1),
     &                      KSTAR(I2), RADIUS(I1), RADIUS(I2), QPERI,
     &                      SEMI, ES0, ECC
   15         FORMAT (' CHAOS => SPIRAL    T KS NAM K* R* QP A E0 E ',
     &                               F9.2,I4,2I6,2I4,1P,4E10.2,0P,2F7.3)
*       Activate spiral indicator and save time, pericentre & eccentricity.
	      KSTAR(I) = -2
	      TOSC(IC) = TIME
	      RP(IC) = QPERI
	      ES(IC) = ECC
              NSP = NSP + 1
              IF (KZ(8).GT.3) THEN
                  CALL BINEV(IPAIR)
              END IF
	      GO TO 30
	  END IF
*
*       Combine the two stars inelastically in case of chaos disruption.
          IF (IDIS.GT.0.AND.QPERI.LT.RCOLL) THEN
              R1 = MAX(RADIUS(I1),RADIUS(I2))
              WHICH1 = ' CHAOS  '
              IF (KSTAR(I).EQ.-2) WHICH1 = ' SPIRAL '
              WRITE (6,20)  WHICH1, IPAIR, NAME(I1), NAME(I2),
     &                      KSTAR(I1), KSTAR(I2), R1, R(IPAIR), QPERI,
     &                      SEMI, ECC, ES0, XN
   20         FORMAT (' DISRUPTED',A8,'  KS NM K* R* R QP A E E0 n ',
     &                                   I4,2I6,2I4,1P,4E10.2,0P,3F7.3)
*       Update chaos variables at end of routine SPIRAL (argument < 0).
              IQCOLL = 1
              IF (KSTAR(I).EQ.-2) IQCOLL = 2
              II = -I
              CALL SPIRAL(II)
              KSTAR(I) = 0
              CALL XVPRED(I,0)
              KSPAIR = IPAIR
              CALL CMBODY(R(IPAIR),2)
              DMR = -1.0
              GO TO 30
          END IF
      ELSE IF (KSTAR(I).EQ.-2) THEN
          I1 = 2*IPAIR - 1
          I2 = I1 + 1
          IF (QPERI.LT.2.0*MAX(RADIUS(I1),RADIUS(I2))) THEN
*       Determine indices for primary & secondary star (donor & accretor).
              J1 = I1
              IF (RADIUS(I2).GT.RADIUS(I1)) J1 = I2
*       Define mass ratio and evaluate Roche radius for the primary.
              Q0 = BODY(J1)/(BODY(I) - BODY(J1))
              Q1 = Q0**0.3333
              Q2 = Q1**2
              RL1 = 0.49*Q2/(0.6*Q2 + LOG(1.0D0 + Q1))*SEMI
*       Check Roche radius but skip RESET call (no second EMERGE correction).
              IF (RADIUS(J1).GT.RL1.AND.IPHASE.NE.7) THEN
                  WRITE (6,25) NAME(I1), NAME(I2), KSTAR(I1), KSTAR(I2),
     &                         ECC, ES0, RCOLL, RL1, QPERI, SEMI, TC
   25             FORMAT (' DISRUPTED SPIRAL    NM K* E E0 RC RL QP A',
     &                    ' TC ',2I6,2I4,2F7.3,1P,5E10.2)
*       Obtain KS variables at pericentre and enforce collision or CE.
                  CALL KSPERI(IPAIR)
                  KSTAR(I) = 0
                  KSPAIR = IPAIR
                  IQCOLL = 2
                  CALL CMBODY(QPERI,2)
                  DMR = -1.0
*       Note that same result achieved by TERM SPIRAL in case IPHASE = 7.
              END IF
          END IF
      END IF
*
   30 RETURN
*
      END
back to top