ksterm.f

Go to the documentation of this file.

      SUBROUTINE ksterm
*
*
*       Termination of KS regularization.
*       ---------------------------------
*
      include 'common6.h'
      common/slow0/  range,islow(10)
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflag(mmax)
      REAL*8  SAVE(15)
*
*
*       Copy pair index from COMMON save and define KS components & c.m.
      ipair = kspair
      i1 = 2*ipair - 1
      i2 = i1 + 1
      icm = n + ipair
      jmin = 0
*
*       Prepare termination at block time (KS, triple, quad, chain or merge).
      IF (time.LE.tblock.AND.iphase.LT.9) THEN
          time0 = tblock
*       Skip KS integration for unperturbed orbit or T0(I1) at end of block.
          IF (list(1,i1).EQ.0.OR.t0(i1).EQ.time0) go to 3
*
*       See whether the time interval should be modified by KSLOW procedure.
          IF (kslow(ipair).GT.1) THEN
              imod = kslow(ipair)
              zmod = float(islow(imod))
          ELSE
              zmod = 1.0
          END IF
*
    1     dt0 = time0 - t0(i1)
*       Integrate up to current block time in case interval is too large.
          IF (dt0.GT.step(i1)) THEN
              time = t0(i1) + step(i1)
              h0(ipair) = h(ipair)
              z = -0.5d0*h(ipair)*dtau(ipair)**2
              CALL stumpf(ipair,z)
              CALL ksint(i1)
              dtu = dtau(ipair)
              step(i1) = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu
     &                                                   + r(ipair))*dtu
              step(i1) = zmod*step(i1)
*       Restrict increase of R for superfast particles in one block-step.
              IF (h(ipair).GT.100.0.AND.r(ipair).GT.rmin) go to 3
              go to 1
          END IF
*
*       Determine the last regularized step by Newton-Raphson iteration.
          dtu = dt0/(r(ipair)*zmod)
          dtu = min(dtu,dtau(ipair))
*       Include rare case of zero interval due to subtracting large values.
          dtu = max(dtu,1.0d-10)
          iter = 0
    2     y0 = dt0 - zmod*((one6*tdot3(ipair)*dtu +
     &                             0.5*tdot2(ipair))*dtu + r(ipair))*dtu
          ypr = -((0.5*tdot3(ipair)*dtu + tdot2(ipair))*dtu + r(ipair))
          ypr = zmod*ypr
          dtu = dtu - y0/ypr
          dt1 = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu +
     &                                                     r(ipair))*dtu
          dt1 = zmod*dt1
          iter = iter + 1
          IF (abs(dt0 - dt1).GT.1.0e-10*step(i1).AND.iter.LT.10) go to 2
*
*       Advance the KS solution to next block time and terminate at TIME0.
          dtau(ipair) = dtu
          step(i1) = dt1
          time = t0(i1) + dt1
          h0(ipair) = h(ipair)
          z = -0.5d0*h(ipair)*dtu**2
          CALL stumpf(ipair,z)
          CALL ksint(i1)
    3     time = time0
*
*       Predict X & XDOT for body #JCOMP (note TIME = TBLOCK if second call).
          IF (jcomp.GE.ifirst) THEN
              CALL xvpred(jcomp,-1)
              IF (gamma(ipair).GT.0.2.AND.jcomp.LE.n) THEN
                  jmin = jcomp
*       Initialize T0, X0 & X0DOT for XVPRED & FPOLY on large perturbation.
                  t0(jcomp) = time
                  DO 4 k = 1,3
                      x0(k,jcomp) = x(k,jcomp)
                      x0dot(k,jcomp) = xdot(k,jcomp)
    4             CONTINUE
              END IF
          END IF
      END IF
*
*       Predict coordinates and evaluate potential energy w.r.t. perturbers.
      CALL ksres(ipair,j1,j2,0.0d0)
      np = list(1,i1)
      DO 5 l = 1,np
          jpert(l) = list(l+1,i1)
    5 CONTINUE
      jlist(1) = i1
      jlist(2) = i2
      CALL nbpot(2,np,pot1)
*
*       Rectify the orbit to yield U & UDOT consistent with binding energy.
      CALL ksrect(ipair)
*
*       Retain final KS variables for explicit restart at merge termination.
      IF (time.LE.tblock.AND.iphase.EQ.6) THEN
          hm(nmerge) = h(ipair)
          DO 6 k = 1,4
              um(k,nmerge) = u(k,ipair)
              umdot(k,nmerge) = udot(k,ipair)
    6     CONTINUE
      END IF
*
*       Check optional diagnostic output for disrupted new hard binary.
      IF (kz(8).EQ.0) go to 10
      IF (list(2,i1+1).NE.0.OR.h(ipair).GT.0.0) go to 10
      IF (gamma(ipair).GT.0.5.AND.jcomp.GT.0.OR.iphase.EQ.7) THEN
          IF (jcomp.EQ.0.OR.iphase.EQ.7) jcomp = i1
          k = 0
          IF (jcomp.GT.n) THEN
              j2 = 2*(jcomp - n)
              k = list(2,j2)
          END IF
          semi = -0.5*body(icm)/h(ipair)
          eb = -0.5*body(i1)*body(i2)/semi
          ri = sqrt((x(1,icm) - rdens(1))**2 +
     &              (x(2,icm) - rdens(2))**2 +
     &              (x(3,icm) - rdens(3))**2)
          WRITE (8,8)  time+toff, name(i1), name(i2), k, name(jcomp),
     &                 body(jcomp), eb, semi, r(ipair), gamma(ipair), ri
    8     FORMAT (' END BINARY   T =',f8.1,'  NAME = ',2i6,i3,i6,
     &                      '  M(J) =',f8.4,'  EB =',f10.5,'  A =',f8.5,
     &                          '  R =',f8.5,'  G =',f5.2,'  RI =',f5.2)
      END IF
*
   10 IF (kz(10).GT.1) THEN
          ri = sqrt((x(1,icm) - rdens(1))**2 +
     &              (x(2,icm) - rdens(2))**2 +
     &              (x(3,icm) - rdens(3))**2)
*       Include error of the bilinear relation (cf. Book Eqn.(4.30)).
          err = u(4,ipair)*udot(1,ipair) - u(3,ipair)*udot(2,ipair) +
     &          u(2,ipair)*udot(3,ipair) - u(1,ipair)*udot(4,ipair)
          WRITE (6,15)  time+toff, body(i1), body(i1+1), dtau(ipair),
     &                  r(ipair), ri, h(ipair), ipair, gamma(ipair),
     &                  step(i1), list(1,i1), list(1,icm), err
   15     FORMAT (/,' END KSREG    TIME =',f8.2,1p,2e9.1,0p,f6.3,1p,
     &                e10.1,0p,f7.2,f9.2,i5,f8.3,1p,e10.1,2i5,1p,e10.2)
      END IF
*
*       Obtain global coordinates & velocities.
      CALL resolv(ipair,2)
*
*       Correct for differential potential energy due to rectification.
      CALL nbpot(2,np,pot2)
*       Add correction term with opposite sign for conservation.
*     ECOLL = ECOLL + (POT2 - POT1)
*     IF (ABS(POT1-POT2).GT.0.0001) WRITE (6,16)  POT1,BE(3),POT1-POT2
*  16 FORMAT (' CORRECT:    POT1 BE3 POT1-POT2  ',2F10.6,F10.6)
*
*       Modify c.m. neighbour radius by density contrast and set new values.
      nnb = list(1,icm) + 1
*       Check predicted neighbour number and form volume ratio.
      nbp = min(alpha*sqrt(float(nnb)*rs(icm))/(rs(icm)**2),znbmax)
      nbp = max(nbp,int(znbmin))
      a0 = float(nbp)/float(nnb)
*       Re-determine neighbour list on zero membership for distant binary.
      IF (list(1,icm).EQ.0) THEN
          rs0 = 0.1*(abs(x(1,icm)) + abs(x(2,icm)))
          CALL nblist(icm,rs0)
          nnb = list(1,icm) + 1
      END IF
*       Copy all neighbours in the case of merger.
      IF (iphase.EQ.6) THEN
          a0 = 1.0
          nbp = nnb - 1
      END IF
      IF (rs(icm).GT.-100.0*body(icm)/h(ipair)) a0 = 1.0
*       Accept old c.m. values for small length scale ratio or H > 0.
      rs(i1) = rs(icm)*a0**0.3333
      rs(i1+1) = rs(i1)
      rs2 = rs(i1)**2
*
*       Select neighbours for components inside the modified c.m. sphere.
   20 nnb1 = 1
      DO 25 l = 2,nnb
          j = list(l,icm)
          rij2 = (x(1,icm) - x(1,j))**2 + (x(2,icm) - x(2,j))**2 +
     &                                    (x(3,icm) - x(3,j))**2
*       Ensure that at least the predicted neighbour number is reached.
          IF (rij2.LT.rs2.OR.(nnb1 + nnb - l.LE.nbp.AND.
     &        nnb1.LT.nnbmax-1)) THEN
              nnb1 = nnb1 + 1
              ilist(nnb1) = j
          END IF
   25 CONTINUE
*
*       Check that there is space for adding dominant component later.
      IF (nnb1.GE.nnbmax.AND.iphase.NE.6) THEN
          rs2 = 0.9*rs2
          go to 20
      END IF
*
*       Reduce pair index, total number & single particle index.
      npairs = npairs - 1
      ntot = n + npairs
      ifirst = 2*npairs + 1
*
*       Save name of components & flag for modifying LISTD in UPDATE.
      jlist(1) = name(i1)
      jlist(2) = name(i1+1)
      jlist(3) = list(2,i1+1)
*
*       Skip adjustment of tables if last or only pair being treated.
      IF (ipair.EQ.npairs + 1) go to 60
*
*       Move the second component before the first.
      DO 50 kcomp = 2,1,-1
          i = 2*ipair - 2 + kcomp
*
          DO 30 k = 1,3
              save(k) = x(k,i)
              save(k+3) = x0dot(k,i)
   30     CONTINUE
*       Current velocity has been set in routine RESOLV.
          save(7) = body(i)
          save(8) = rs(i)
          save(9) = radius(i)
          save(10) = tev(i)
          save(11) = body0(i)
          save(12) = tev0(i)
          save(13) = epoch(i)
          save(14) = spin(i)
          save(15) = zlmsty(i)
          namei = name(i)
          ksi = kstar(i)
          last = 2*npairs - 1 + kcomp
*
*       Move up global variables of other components.
          DO 40 j = i,last
              j1 = j + 1
              DO 35 k = 1,3
                  x(k,j) = x(k,j1)
*       Copy latest X & X0DOT (= 0) of single components for predictor.
                  x0(k,j) = x(k,j)
                  x0dot(k,j) = x0dot(k,j1)
   35         CONTINUE
              body(j) = body(j1)
              rs(j) = rs(j1)
              radius(j) = radius(j1)
              tev(j) = tev(j1)
              tev0(j) = tev0(j1)
              body0(j) = body0(j1)
              epoch(j) = epoch(j1)
              spin(j) = spin(j1)
              zlmsty(j) = zlmsty(j1)
              name(j) = name(j1)
              kstar(j) = kstar(j1)
              step(j) = step(j1)
              t0(j) = t0(j1)
              k = list(1,j1) + 1
              IF (k.EQ.1) k = 2
*       Transfer unmodified neighbour lists (include flag in 2nd comp).
              DO 38 l = 1,k
                  list(l,j) = list(l,j1)
   38         CONTINUE
   40     CONTINUE
*
*       Set new component index and copy basic variables.
          i = last + 1
          DO 45 k = 1,3
              x(k,i) = save(k)
              x0dot(k,i) = save(k+3)
              xdot(k,i) = save(k+3)
   45     CONTINUE
          body(i) = save(7)
          rs(i) = save(8)
          radius(i) = save(9)
          tev(i) = save(10)
          body0(i) = save(11)
          tev0(i) = save(12)
          epoch(i) = save(13)
          spin(i) = save(14)
          zlmsty(i) = save(15)
          name(i) = namei
          kstar(i) = ksi
   50 CONTINUE
*
*       Include removal of the circularization name NAMEC from chaos table.
      IF (kstar(icm).GE.10.AND.nchaos.GT.0.AND.iphase.NE.6) THEN
*       Note that NAMEC may remain dormant for hierarchical systems.
          ii = -icm
          CALL spiral(ii)
      END IF
*
*       Update all regularized variables.
      CALL remove(ipair,2)
*
*       Remove old c.m. from all COMMON tables (no F & FDOT correction).
      CALL remove(icm,3)
*
*       Set new global index of first & second component.
   60 icomp = 2*npairs + 1
      jcomp = icomp + 1
*
*       Save c.m. neighbour list for routine FPOLY1/2 (may be renamed below).
      ilist(1) = nnb1 - 1
      DO 65 l = 1,nnb1
          list(l,icomp) = ilist(l)
   65 CONTINUE
*
*       Modify all relevant COMMON list arrays.
      CALL update(ipair)
*
*       Check replacing of single KS component by corresponding c.m.
   70 IF (list(2,icomp).LT.icomp) THEN
          j2 = list(2,icomp)
          j = kvec(j2) + n
          DO 80 l = 2,nnb1
              IF (l.LT.nnb1.AND.list(l+1,icomp).LT.j) THEN
                  list(l,icomp) = list(l+1,icomp)
              ELSE
                  list(l,icomp) = j
              END IF
   80     CONTINUE
*       Check again until first neighbour > ICOMP.
          go to 70
      END IF
*
*       Make space for dominant component and copy members to JCOMP list.
      DO 90 l = nnb1,2,-1
          list(l+1,icomp) = list(l,icomp)
          list(l+1,jcomp) = list(l,icomp)
   90 CONTINUE
*
*       Set dominant component in first location and specify membership.
      list(2,icomp) = jcomp
      list(2,jcomp) = icomp
      list(1,icomp) = nnb1
      list(1,jcomp) = nnb1
*
*       Initialize T0, T0R, X0 & X0DOT for both components.
      t0(icomp) = time
      t0(jcomp) = time
      t0r(icomp) = time
      t0r(jcomp) = time
      DO 95 k = 1,3
          x0(k,icomp) = x(k,icomp)
          x0(k,jcomp) = x(k,jcomp)
          x0dot(k,icomp) = xdot(k,icomp)
          x0dot(k,jcomp) = xdot(k,jcomp)
   95 CONTINUE
*
*       Form new force polynomials (skip triple, quad, merge & chain).
      IF (iphase.LT.4) THEN
*       Predict current coordinates & velocities for the neighbours.
          CALL xvpred(icomp,nnb1)
*
*       Obtain new polynomials & steps.
          CALL fpoly1(icomp,jcomp,2)
          CALL fpoly2(icomp,jcomp,2)
*
*       Improve force polynomials of strong perturber after rectification.
          IF (jmin.GE.ifirst) THEN
              CALL fpoly1(jmin,jmin,0)
              CALL fpoly2(jmin,jmin,0)
          END IF
      END IF
*
*       Check updating of global index for chain c.m.
      IF (nch.GT.0) THEN
          CALL chfind
      END IF
*
      RETURN
*
      END