recoil.f

Go to the documentation of this file.

      SUBROUTINE recoil(IEND)
*
*
*       Binary analysis and final output.
*       ---------------------------------
*
      include 'commonc.h'
      include 'common2.h'
      REAL*8  r2(nmx,nmx),rc(3),vc(3),xx(3,3),vv(3,3)
      INTEGER  ij(nmx)
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      common/cpert/  rgrav,gpert,ipert,npert
      common/calls/  tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
      common/ccoll2/  qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
     &                ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
      common/clump/   bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(nmx,5),isys(5)
      SAVE
*
*
*       Sort particle separations (I1 & I2 form closest pair).
      CALL r2sort(ij,r2)
      i1 = ij(1)
      i2 = ij(2)
      i3 = ij(3)
      i4 = ij(4)
      IF (n.EQ.2) THEN
          i3 = i2
          i4 = i1
      ELSE IF (n.EQ.3) THEN
          i4 = i1
      ELSE IF (n.GE.4) THEN
*       Determine indices of second closest pair (avoid pair I1-I2).
          rx1 = 1.0
          rx0 = r2(i1,i2)
          DO 2 j1 = 1,n
              IF (j1.EQ.i1.OR.j1.EQ.i2) go to 2
              DO 1 j2 = j1+1,n
                  IF (j2.EQ.i1.OR.j2.EQ.i2) go to 1
                  IF (r2(j1,j2).LT.rx1.AND.r2(j1,j2).GT.rx0) THEN
                      rx1 = r2(j1,j2)
                      i3 = j1
                      i4 = j2
                  END IF
    1         CONTINUE
    2     CONTINUE
      END IF
      k1 = i3
      k2 = i4
*
*       Ensure that original KS pair is considered (H > 0 is possible).
      IF (n.LE.3.AND.iend.EQ.0) THEN
          DO 5 k = 1,n
              IF (iname(k).EQ.1) k1 = k
              IF (iname(k).EQ.2) k2 = k
    5     CONTINUE
      END IF
*
*       Form output diagnostics.
      vrel2 = 0.0d0
      vrel21 = 0.0d0
      vrel23 = 0.0d0
      rdot = 0.0d0
      rc2 = 0.0d0
      vc2 = 0.0d0
      rcdot = 0.0d0
      mb = m(i1) + m(i2)
*
      DO 10 k = 1,3
          j1 = 3*(i1-1) + k
          j2 = 3*(i2-1) + k
          j3 = 3*(i3-1) + k
          j4 = 3*(i4-1) + k
          l1 = 3*(k1-1) + k
          l2 = 3*(k2-1) + k
          vrel2 = vrel2 + (v(j1) - v(j2))**2
          vrel21 = vrel21 + (v(j3) - v(j4))**2
          vrel23 = vrel23 + (v(l1) - v(l2))**2
          rdot = rdot + (x(j1) - x(j2))*(v(j1) - v(j2))
          rc(k) = (m(i1)*x(j1) + m(i2)*x(j2))/mb
          vc(k) = (m(i1)*v(j1) + m(i2)*v(j2))/mb
          rc2 = rc2 + (rc(k) - x(j3))**2
          vc2 = vc2 + (vc(k) - v(j3))**2
          rcdot = rcdot + (rc(k) - x(j3))*(vc(k) - v(j3))
          xx(k,1) = x(j1)
          xx(k,2) = x(j2)
          xx(k,3) = x(j3)
          vv(k,1) = v(j1)
          vv(k,2) = v(j2)
          vv(k,3) = v(j3)
   10 CONTINUE
*
*       Evaluate orbital elements.
      rb = sqrt(r2(i1,i2))
      rb1 = sqrt(r2(i3,i4))
      r13 = sqrt(r2(i1,i3))
      r24 = sqrt(r2(i2,i4))
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
      ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
      eb = -m(i1)*m(i2)/(2.0d0*semi)
      et = -m(i1)*m(i2)/(2.0d0*semi*cm(8))
*
*       Determine parameters for possible hierarchy of body #I3.
      rcp = sqrt(rc2)
      a1 = 2.0/rcp - vc2/(mb + m(i3))
      a1 = 1.0/a1
      ecc1 = sqrt((1.0 - rcp/a1)**2 + rcdot**2/(a1*(mb + m(i3))))
      pmin = a1*(1.0 - ecc1)
*
*       Obtain binding energy of K1 & K2 (in case I1 & I2 not bound).
      semi2 = 2.0d0/sqrt(r2(k1,k2)) - vrel23/(m(k1) + m(k2))
      semi2 = 1.0/semi2
      eb2 = -m(k1)*m(k2)/(2.0d0*semi2)
*
*       Include possible second binary (suppress positive energy).
      IF (n.GT.3) THEN
          mb1 = m(i3) + m(i4)
          semi1 = 2.0d0/rb1 - vrel21/mb1
          semi1 = 1.0/semi1
          eb1 = -m(i3)*m(i4)/(2.0d0*semi1)
          eb1 = min(eb1,0.0d0)
      ELSE
          r24 = 1.0e+10
          IF (timec.LE.0.0d0) eb1 = 0.0
      END IF
*
*       Ensure that perturbed boundary exceeds system size.
      rmax = max(sqrt(r2(i1,i3)),sqrt(r2(i2,i4)))
      rmaxc = max(rmaxc,1.2*rmax)
*
*       Save initial energies (note: K1 & K2 may be most energetic).
      IF (iend.EQ.0) THEN
          e0 = energy
          de = 0.0
*         IF (TIMEC.GT.0.0D0) GO TO 50
          ebch0 = cm(9)
          eb0 = min(eb,eb2)
          eb10 = eb1
          a0 = semi
          ecc0 = ecc
          IF (eb2.LT.eb) a0 = semi2
*       Note closest pair (#I1/I2 may be intruder or wide binary at peri).
          IF (n.EQ.3) THEN
              name1 = namec(1)
              name2 = namec(2)
          ELSE
              name1 = namec(i1)
              name2 = namec(i2)
          END IF
          go to 50
      ELSE IF (iend.EQ.1) THEN
          de = de + abs(energy - e0)
          e0 = energy
*       Initialize second binary on subsequent absorption.
          IF (eb10.EQ.0.0d0) eb10 = eb1
          go to 50
      END IF
*
*       Determine nearest and most distant binary perturber (N = 4).
      IF (r13.LT.r24) THEN
          im = i3
          rm = r13
          imax = i4
          rmax = r24
      ELSE
          im = i4
          rm = r24
          imax = i3
          rmax = r13
      END IF
*
*       Estimate perturbations on smallest binary.
      IF (n.GT.3) THEN
          IF (n.GT.4) rm = max(r13,r24)
          gb = 2.0*mb1/mb*(rb/rm)**3
          IF (eb1.LT.0.0) gb = gb*m(im)/mb1
          g4 = 2.0*m(imax)/(mb + m(im))*(rm/rmax)**3
      ELSE
          gb = m(i3)/mb*(rb/r13)**3
          g4 = 0.0
          eb1 = 0.0
      END IF
*
*       Form net binary energy change (added to CHCOLL in routine CHTERM).
      IF (iend.EQ.2) THEN
          cm(9) = eb + eb1 - cm(9) + ecoll1
      END IF
*
*       Set relative energy production and total energy change.
      db = cm(9)/ebch0
      de = de + abs(energy - e0)
*
*       Provide diagnostics for exchange (membership may have changed).
      IF (name1 + name2.NE.namec(i1) + namec(i2).AND.eb.LT.0.0) THEN
          isub = isys(5)
          tch = t0s(isub) + timec
          WRITE (6,15)  tch, name1, name2, namec(i1), namec(i2), ecc0,
     &                  ecc, a0, semi, eb0, eb
   15     FORMAT (' EXCHANGE    T NAM E0 E A0 A EB0 EB ',
     &                          f9.2,4i6,2f7.3,1p,2e9.1,2e10.1)
      END IF
*
*       Print final binary if relative energy increase exceeds 0.1.
      IF (db.GT.0.1.AND.semi.GT.0.0) THEN 
*       Scale binding energy of relative motion by initial value.
          e1 = (energy - eb - eb1)/eb0
          eb = eb/energy
          eb1 = eb1/energy
          WRITE (6,20)  namec(i1), namec(i2), semi, ecc, eb, gb, g4,
     &                  eb1, e1, et, db
   20     FORMAT (' CHAIN RECOIL','  NAM =',2i6,'  A =',1p,e8.1,
     &            '  E =',0p,f5.2,'  EB =',f5.2,'  GB =',1p,e8.1,
     &            '  G4 =',e8.1,'  EB1 =',0p,f5.2,'  E1 =',f5.2,
     &            '  ET =',f6.3,'  DB =',f5.1)
      END IF
*
*       Include output of strong interaction from chaos case.
      IF (iend.EQ.2.AND.pmin.LT.2.0*semi.AND.db.GT.0.1) THEN
          CALL inclin(xx,vv,rc,vc,alpha)
          WRITE (6,25)  namec(i3), ecc, ecc1, pmin/semi, rcdot/rcp,
     &                  semi, a1, rcp, gb, 180.0*alpha/3.14
   25     FORMAT (' RECOIL:    NAM E E1 PM/A RD A A1 RP GB IN ',
     &                         i6,f8.4,f6.2,f5.1,f6.1,1p,4e10.2,0p,f7.1)
      END IF
*
*       Check optional print diagnostics of chain regularization.
      IF (kz30.GT.1.AND.iend.EQ.2) THEN
          tcr = mass**2.5/abs(2.0d0*energy)**1.5
          tc = timec/tcr
          ec = energy/cm(8)
          WRITE (6,30)  i1, i2, i3, i4, rb, r13, r24, de, tc, nstep1,
     &                  nreg, npert, db, ec
   30     FORMAT (/,' END CHAIN  ',4i3,'  RB =',1pe8.1,'  R13 =',e8.1,
     &              '  R24 =',e8.1,'  DE =',e8.1,'  TC =',0p,f5.1,'  #',
     &                 i5,i4,i3,'  DB =',f5.2,'  EC =',f6.3)
      END IF
*
   50 RETURN
*
      END