ksinit.f

Go to the documentation of this file.

      SUBROUTINE ksinit
*
*
*       Initialization of KS regularization.
*       ------------------------------------
*
      include 'common6.h'
      REAL*8  q(3),rdot(3),ui(4),vi(4),a1(3,4)
      INTEGER ipipe(9)
      SAVE  ipipe
      DATA  ipipe /9*0/
*
*
*       Set new global indices of the components and current pair index.
      icomp = 2*npairs - 1
      jcomp = icomp + 1
      ipair = npairs
*
*       Add body #N in case the only neighbour was removed in KSREG.
      IF (list(1,icomp).EQ.0) THEN
          list(2,icomp) = n
          list(2,jcomp) = n
          list(1,icomp) = 1
          list(1,jcomp) = 1
      END IF
*
*       Specify mass, neighbour radius, name, radius & type for new c.m.
      body(ntot) = body(icomp) + body(jcomp)
      rs(ntot) = rs(icomp)
      name(ntot) = nzero + name(icomp)
      radius(ntot) = 0.0
      tev(ntot) = 1.0e+10
      tev0(ntot) = 1.0e+10
      body0(ntot) = body(ntot)
      epoch(ntot) = time*tstar
      kstar(ntot) = 0
      imod = 1
*
*       Define c.m. coordinates & velocities and set XDOT for components.
      DO 10 k = 1,3
          x(k,ntot) = (body(icomp)*x(k,icomp) + body(jcomp)*x(k,jcomp))/
     &                                                        body(ntot)
          x0dot(k,ntot) = (body(icomp)*x0dot(k,icomp) + body(jcomp)*
     &                                        x0dot(k,jcomp))/body(ntot)
          xdot(k,ntot) = x0dot(k,ntot)
          xdot(k,icomp) = x0dot(k,icomp)
          xdot(k,jcomp) = x0dot(k,jcomp)
   10 CONTINUE
*
*       Obtain force polynomial for c.m. with components ICOMP & JCOMP.
      nnb = list(1,icomp)
*
*       Predict current coordinates & velocities for the neighbours.
      CALL xvpred(icomp,nnb)
*
*       Obtain new polynomials & steps (first F & FDOT, then F2DOT & F3DOT).
      CALL fpoly1(icomp,jcomp,1)
      CALL fpoly2(ntot,ntot,1)
*
*       Skip KS initialization at merger termination (H, U & UDOT in RESET).
      IF (iphase.EQ.7) THEN
          eb = 2.0*ebh
          go to 50
      END IF
*
*       Define relative coordinates and velocities in physical scaled units.
      DO 20 k = 1,3
          q(k) = x(k,icomp) - x(k,jcomp)
          rdot(k) = x0dot(k,icomp) - x0dot(k,jcomp)
   20 CONTINUE
*
*       Introduce regularized variables using definition of Book.
      r(ipair) = sqrt(q(1)**2 + q(2)**2 + q(3)**2)
*
*       Initialize the regularized coordinates according to sign of Q(1).
      IF (q(1).LE.0.0d0) THEN
          ui(3) = 0.0d0
          ui(2) = sqrt(0.5d0*(r(ipair) - q(1)))
          ui(1) = 0.5d0*q(2)/ui(2)
          ui(4) = 0.5d0*q(3)/ui(2)
      ELSE
          ui(4) = 0.0d0
          ui(1) = sqrt(0.5d0*(r(ipair) + q(1)))
          ui(2) = 0.5d0*q(2)/ui(1)
          ui(3) = 0.5d0*q(3)/ui(1)
      END IF
*
*       Set current transformation matrix.
      CALL matrix(ui,a1)
*
*       Form regularized velocity and set initial KS coordinates.
      tdot2(ipair) = 0.0
      DO 30 k = 1,4
          udot(k,ipair) = 0.50d0*(a1(1,k)*rdot(1) + a1(2,k)*rdot(2) +
     &                                              a1(3,k)*rdot(3))
*       Note that A1(J,K) is the transpose of A1(K,J).
          u(k,ipair) = ui(k)
          u0(k,ipair) = u(k,ipair)
          tdot2(ipair) = tdot2(ipair) + 2.0d0*ui(k)*udot(k,ipair)
   30 CONTINUE
*
*       Evaluate initial binding energy per unit mass (singular form).
      h(ipair) = (2.0d0*(udot(1,ipair)**2 + udot(2,ipair)**2 +
     &                   udot(3,ipair)**2 + udot(4,ipair)**2) -
     &                                              body(ntot))/r(ipair)
      eb = h(ipair)*body(icomp)*body(jcomp)/body(ntot)
*
*       Form perturber list.
   50 CALL kslist(ipair)
*
*       Transform any unperturbed hard binary to apocentre and set time-step.
      semi = -0.5*body(ntot)/h(ipair)
      IF (list(1,icomp).EQ.0.AND.eb.LT.ebh.AND.semi.LT.rmin) THEN
          tk = twopi*semi*sqrt(semi/body(ntot))
*       Note TIME is not commensurate after KSPERI (cf. CHTERM & STEPS).
          IF (iphase.NE.7.AND.iphase.NE.8) THEN
              DO 55 k = 1,4
                  vi(k) = udot(k,ipair)
   55         CONTINUE
*       Determine pericentre time (TP < 0 if TDOT2 < 0) and add TK/2.
              CALL tperi(semi,ui,vi,body(ntot),tp)
*       Note: apocentre to apocentre gives almost zero step.
              step(icomp) = 0.5*min(tk,step(ntot)) - tp
*       Transform KS variables to peri and by pi/2 to apocentre (skip apo).
              IF (abs(tdot2(ipair)).GT.1.0e-12.OR.r(ipair).LT.semi) THEN
                  time0 = time
                  CALL ksperi(ipair)
                  CALL ksapo(ipair)
*       Reset TIME to quantized value (small > 0 or < 0 possible initially).
                  time = time0
                  time = max(time,0.0d0)
              ELSE IF (tdot2(ipair).GT.0.0) THEN
                  tdot2(ipair) = -1.0e-20
              END IF
          END IF
      END IF
*
*       Estimate an appropriate KS slow-down index for G < GMIN.
      IF (list(1,icomp).EQ.0.AND.semi.GT.0.0) THEN
          tk = twopi*semi*sqrt(semi/body(ntot))
          IF (kz(26).GT.0.AND.step(ntot).GT.tk) THEN
              imod = 1 + log(step(ntot)/tk)/0.69
              imod = min(imod,5)
          END IF
      END IF
*
*       Specify zero membership and set large steps for second component.
      list(1,jcomp) = 0
*       Set large step for second component to avoid detection.
      step(jcomp) = 1.0e+06
      stepr(jcomp) = 1.0e+06
*
*       Obtain polynomials for perturbed KS motion (standard case & merger).
      CALL kspoly(ipair,imod)
*
*       Obtain apocentre distance.
      semi = -0.5*body(ntot)/h(ipair)
      ecc2 = (1.0-r(ipair)/semi)**2 + tdot2(ipair)**2/(body(ntot)*semi)
      rap = semi*(1.0 + sqrt(ecc2))
*
*       Include suggestion for monitoring hyperbolic encounters (suppressed).
*     IF (SEMI.LT.0.0) THEN
*         PMIN = SEMI*(1.0 - SQRT(ECC2))
*         WRITE (31,56)  TIME+TOFF, NAME(ICOMP), NAME(JCOMP), PMIN
*  56     FORMAT (' HYPERBOLIC    T NAM PMIN ',F8.2,2I6,1P,E10.2)
*     END IF
*
*       Set 2*SEMI as termination scale for hard binary if 2*SEMI < RS/20.
      IF (eb.LT.ebh.AND.rap.LT.0.05*rs(ntot)) THEN 
          r0(ipair) = max(rmin,-body(ntot)/h(ipair))
          r0(ipair) = min(r0(ipair),5.0*rmin)
      ELSE
          r0(ipair) = r(ipair)
      END IF
*
*       Increase regularization counters (#9 & NKSHYP for hyperbolic orbits).
      nksreg = nksreg + 1
      IF (h(ipair).GT.0.0) nkshyp = nkshyp + 1
*
*       Check optional output for new KS.
      IF (kz(10).GT.0) THEN
          ri = sqrt((x(1,ntot) - rdens(1))**2 +
     &              (x(2,ntot) - rdens(2))**2 +
     &              (x(3,ntot) - rdens(3))**2)
          WRITE (6,60)  time+toff, name(icomp), name(jcomp),dtau(ipair),
     &                  r(ipair), ri, h(ipair), ipair, gamma(ipair),
     &                  step(ntot), list(1,icomp), list(1,ntot)
   60     FORMAT (/,' NEW KSREG    TIME =',f8.2,2i6,f12.3,1p,e10.1,
     &                                0p,f7.2,f9.2,i5,f8.3,1p,e10.1,2i5)
      END IF
*
*       Include limited diagnostics for NS or BH hard binary formation.
      IF (max(kstar(icomp),kstar(jcomp)).GE.13.AND.eb.LT.ebh.AND.
     &    iphase.NE.7) THEN
          id = 0
          np = ipipe(1)
*       See whether at least one component was recorded in previous pair.
          DO 62 j = 2,np+1
              IF (ipipe(j).EQ.name(icomp).OR.
     &            ipipe(j).EQ.name(jcomp)) id = id + 1
   62     CONTINUE
*       Print diagnostics if NS/BH binary not listed among four last pairs.
          IF (id.LE.1) THEN
              pd = days*semi*sqrt(abs(semi)/body(ntot))
              WRITE (6,65)  time+toff, name(icomp), name(jcomp),
     &                      kstar(icomp), kstar(jcomp), kstar(ntot),
     &                      sqrt(ecc2), pd, eb
   65         FORMAT (' NS/BH BINARY    T NM K* E P EB ',
     &                                  f8.1,2i6,3i4,f7.3,1p,e9.1,e11.2)
          END IF
*       Update table and membership by removing first two.
          IF (np.GE.8) THEN
*       Note there are at most 4 pairs with entries in IPIPE(2:9).
              DO 66 j = 2,6
                  ipipe(j) = ipipe(j+2)
                  ipipe(j+1) = ipipe(j+3)
   66         CONTINUE
              np = np - 2
          END IF
*       Add NAME of each component in NP+2/3 and increase membership by 2.
          ipipe(np+2) = name(icomp)
          ipipe(np+3) = name(jcomp)
          ipipe(1) = np + 2
      END IF
*
*       Modify the termination criterion according to value of NPAIRS.
      IF (npairs.GT.kmax - 3) gmax = 0.8*gmax
      IF (npairs.LT.kmax - 5.AND.gmax.LT.0.001) gmax = 1.2*gmax
      IF (npairs.EQ.kmax) WRITE (6,70)  npairs, time+toff
   70 FORMAT (5x,'WARNING!   MAXIMUM KS PAIRS   NPAIRS TIME',i5,f9.2)
*
*       Initialize prediction variables to avoid skipping KS components.
      DO 75 kcomp = 1,2
          jdum = 2*npairs - 2 + kcomp
          DO 72 k = 1,3
              x0(k,jdum) = x(k,jdum)
              x0dot(k,jdum) = 0.0d0
              f(k,jdum) = 0.0d0
              fdot(k,jdum) = 0.0d0
              d2(k,jdum) = 0.0d0
              d3(k,jdum) = 0.0d0
              d2r(k,jdum) = 0.0d0
              d3r(k,jdum) = 0.0d0
   72     CONTINUE
   75 CONTINUE
*
*       See whether either component has been regularized recently.
      nnb = listd(1) + 1
      k = 0
*       Check case of initial binary and loop over disrupted pairs.
      IF (iabs(name(icomp) - name(jcomp)).EQ.1) THEN
          IF (name(icomp).LE.2*nbin0) k = -1
      END IF
      DO 80 l = 2,nnb
          IF (name(icomp).EQ.listd(l).OR.name(jcomp).EQ.listd(l)) k = -1
   80 CONTINUE
      IF (h(ipair).GT.0.0) k = 0
*
*       Treat mergers as new binaries and ensure chain/hard KS as standard.
      IF (iphase.EQ.6) THEN
          k = 0
      ELSE IF (k.EQ.0) THEN
          IF (iphase.EQ.8.OR.eb.LT.ebh) k = -1
      END IF
*       Set flag to distinguish between existing and new binaries (K = -1/0).
      list(2,jcomp) = k
*
*       Check diagnostic output of new hard binary.
      IF (kz(8).GT.0.AND.k.EQ.0) THEN
          IF (eb.GT.ebh) go to 100
          semi = -0.5*body(ntot)/h(ipair)
          ri = sqrt((x(1,ntot) - rdens(1))**2 +
     &              (x(2,ntot) - rdens(2))**2 +
     &              (x(3,ntot) - rdens(3))**2)
          WRITE (8,90)  time+toff, name(icomp), name(jcomp), k,
     &                  body(icomp),body(jcomp), eb, semi, r(ipair),
     &                  gamma(ipair), ri
   90     FORMAT (' NEW BINARY   T =',f7.1,'  NAME = ',2i6,i3,
     &                        '  M =',1p,2e9.1,'  EB =',e9.1,
     &                        '  A =',e9.1,'  R =',e9.1,'  G =',e9.1,
     &                        '  RI =',e9.1)
          CALL flush(8)
      END IF
*
  100 RETURN
*
      END