binout.f

Go to the documentation of this file.

      SUBROUTINE binout
*
*
*       Binary analysis & output.
*       -------------------------
*
      include 'common6.h'
      common/echain/  ech
      INTEGER  jd(15),jeb(15),je(10)
*
*
*       Define semi-major axis & binding energy of hard binary (= KT).
      a0 = 1.0/float(nzero)
      eb0 = -0.25/float(nzero)
*       Initialize counters & variables.
      DO 10 j = 1,15
          jd(j) = 0
          jeb(j) = 0
   10 CONTINUE
      DO 20 j = 1,10
          je(j) = 0
   20 CONTINUE
      ebmax = 0.0
      disp = 0.0
      emax = 0.0
      jor = 0
      jex = 0
      jc = 0
      jlast = 0
      klast = 0
      e(1) = 0.0d0
      e(2) = 0.0d0
      npop(1) = 0
      npop(2) = 0
      npop(3) = n - 2*npairs
      newb = 0
      k10 = 0
*
*       Obtain relevant distributions of KS binaries.
      DO 50 j = 1,npairs
          i = n + j
          IF (h(j).GE.0.0.OR.body(i).EQ.0.0d0) go to 50
          j1 = 2*j - 1
          j2 = 2*j
*
*       Count original & exchanged binaries and core members.
          IF (iabs(name(j1) - name(j2)).EQ.1) THEN
              IF (name(j1).LE.2*nbin0) jor = jor + 1
          ELSE
              IF (min(name(j1),name(j2)).LE.2*nbin0) jex = jex + 1
          END IF
          ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                   (x(3,i) - rdens(3))**2
          IF (ri2.LT.rc2) jc = jc + 1
*
*       Adopt logarithmic distribution of semi-major axis (factor 2).
          semi = -0.5d0*body(i)/h(j)
          IF (semi.LT.a0) THEN
              k = 2 + log10(a0/semi)/log10(2.0)
              k = min(k,10)
              jlast = max(jlast,k)
              jd(k) = jd(k) + 1
          ELSE
              jd(1) = jd(1) + 1
          END IF
*
*       Form eccentricity dispersion & histogram.
          ecc2 = (1.0 - r(j)/semi)**2 + tdot2(j)**2/(body(i)*semi)
          disp = disp + ecc2
          IF (ecc2.GT.emax) emax = ecc2
          k = 1 + 10.0*sqrt(ecc2)
          IF (k.LE.10) je(k) = je(k) + 1
*
*       Set up logarithmic distribution of binding energy (factor 2).
          eb = -0.5*body(j1)*body(j2)/semi
          IF (eb.LT.eb0) THEN
              k = 2 + log10(eb/eb0)/log10(2.0)
              k = min(k,14)
              klast = max(klast,k)
              jeb(k) = jeb(k) + 1
              ebmax = min(ebmax,eb)
          ELSE
              jeb(1) = jeb(1) + 1
          END IF
*
*       Define flag to distinguish primordial or non-primordial binary.
          ip = 1
          IF (list(2,j2).EQ.0) ip = 2
          IF (ip.EQ.2) newb = newb + 1
*
*       Sum the respective energies & populations.
          e(ip) = e(ip) + eb
          npop(ip) = npop(ip) + 1
*
*       Perform consistency check in case of over-writing.
          IF (list(2,j2).NE.-1.AND.list(2,j2).NE.0) THEN
              WRITE (6,35)  j, list(1,j1), list(1,j2), list(2,j2),
     &                      name(n+j)
   35         FORMAT (/,5x,'WARNING!   FLAG PROBLEM   PAIR NP NP2 FLAG',
     &                                                     ' NAME ',5i6)
          END IF
*
*       Produce special diagnostics for new binaries (unit 18).
          IF (ip.EQ.2.AND.h(j).LT.0.0) THEN
              vr = ((x(1,i) - rdens(1))*xdot(1,i) +
     &              (x(2,i) - rdens(2))*xdot(2,i) +
     &              (x(3,i) - rdens(3))*xdot(3,i))/sqrt(ri2)
              k = 0
              IF (name(i).LT.0) k = -1
              gx = gamma(j)*(semi*(1.0 + sqrt(ecc2))/r(j))**3
              WRITE (18,40)  ttot, name(j1), name(j2), list(2,j2), k,
     &                       body(j1), body(j2), eb, semi, sqrt(ecc2),
     &                       gx, sqrt(ri2), vr
   40         FORMAT (' T =',f7.1,'  NAME = ',2i6,2i3,'  M =',2f9.4,
     &                       '  EB =',f10.5,'  A =',f8.5,'  E =',f5.2,
     &                      '  GX =',f6.3,'  RI =',f6.2,'  VR =',f4.1)
              CALL flush(18)
          END IF
*       Reset c.m. Roche flag to standard type for non-zero eccentricity.
          IF (kstar(i).GE.10.AND.ecc2.GT.4.0d-06) kstar(i) = 0
          IF (kstar(i).GE.10) k10 = k10 + 1
   50 CONTINUE
*
*       Set fractional gain of binding energy (initial energy = -0.25).
      IF (time.LE.0.0d0) THEN
          jor = nbin0
          db = 0.0
      ELSE
          e(9) = emerge
          db = -4.0*(e(1) + e(2) + e(5) + e(7) + e(9) + ecoll + emdot
     &                                                + ekick - ebin0)
      END IF
*
*       Print eccentricity & energy distributions.
      IF (npairs.EQ.0) go to 70 
      WRITE (6,60)  jor, jex, db, sbcoll, bbcoll, chcoll, jc, nchaos,
     &              newb, k10, (jd(j),j=1,jlast)
   60 FORMAT (/,' BINARIES',4x,'OR =',i5,'  EX =',i3,'  DB =',f7.3,
     &          '  SB =',f8.4,'  BB =',f8.4,'  CH =',f8.4,'  NC =',i3,
     &         '  NCH =',i4,'  NEWB =',i4,'  CIRC =',i4,'  N(A) =',10i4)
*
      IF (disp.GT.0.0d0) disp = sqrt(disp/float(npairs))
      emax = sqrt(emax)
      WRITE (6,65)  disp, emax, (npop(j),j=1,8), (jeb(k),k=1,klast)
   65 FORMAT (' <E> =',f5.2,'  EMAX =',f7.4,'  NPOP =',i5,i5,2i6,i4,3i3,
     &                                                 '  EB/KT =',14i4)
*
*       Form the basic internal energy (binaries & single particles).
   70 etot = 0.0d0
      e(3) = zkin - pot + etide
      DO 80 j = 1,3
          etot = etot + e(j)
   80 CONTINUE
*
*       Sum all separate energy components to monitor conservation.
      etot = etot + esub + emerge + emdot + ecdot + ecoll
*
*       Include energy of chain if NCH > 0.
      IF (nch.GT.0) THEN
          etot = etot + ech
      END IF
*
*       Form the net energy gain in binary interactions (also in events.f).
      degrav = ebin + esub + ebesc + emesc + emerge + egrav - ebin0
      WRITE (6,90)  (e(j),j=1,10), etot, detot, degrav
   90 FORMAT (' ENERGIES   ',10f10.5,'  ETOT =',f11.6,'  DETOT =',f10.6,
     &                                                '  DEGRAV =',f7.3)
*
      RETURN
*
      END