binev.f

Go to the documentation of this file.

      SUBROUTINE binev(IPAIR)
*
*
*       Binary evolution.
*       -----------------
*
      include 'common6.h'
      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),iflagm(mmax)
      REAL*8  m1,m2
      LOGICAL  first
      SAVE  first,iprev,iskip
      DATA  first,iprev,iskip  /.true.,0,0/
*
*
*       Open unit #17 the first time.
      IF (first) THEN
          OPEN (unit=17,status='NEW',form='FORMATTED',file='BINEV')
          first = .false.
          WRITE (17,5)
    5     FORMAT ('   TPHYS  NAM1  NAM2  K1  K2  KC   M1   M2     R1',
     &            '    R2   RI    ECC    SEMI   PERIOD  C',/)
      END IF
*
*       Define KS & c.m. indices and set component masses & radii.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      i = n + ipair
      m1 = body(i1)*smu
      m2 = body(i2)*smu
      r1 = radius(i1)*su
      r2 = radius(i2)*su
*
*       Skip identical or alternating sequences (i.e. KSTAR(I) = -1 & 0).
      isum = kstar(i1) + kstar(i2) + kstar(i)
      IF (isum.EQ.iprev.OR.isum.EQ.iskip) go to 20
      iskip = iprev
      iprev = isum
*
*       Form basic two-body elements (distinguish mergers & ghosts).
      IF (name(i).GT.0.AND.body(i).GT.0.0d0) THEN
          semi = -0.5*body(i)/h(ipair)
          bodyc = body(i)
          ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(bodyc*semi)
      ELSE
          CALL findj(i1,j,im)
          IF (j.LE.0) RETURN
          bodyc = cm(1,im) + cm(2,im)
          semi = -0.5*bodyc/hm(im)
          rj = um(1,im)**2 + um(2,im)**2 + um(3,im)**2 + um(4,im)**2
          ecc2 = (1.0 - rj/semi)**2
      END IF
*
      ecc = sqrt(ecc2)
      p = days*semi*sqrt(abs(semi)/bodyc)
      p = min(p,99999.9d0)
      p = max(p,-1.0d0)
      a0 = min(semi*su,9999.9d0)
      a0 = max(a0,-1.0d0)
      ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                               (x(3,i) - rdens(3))**2
      ri = min(sqrt(ri2),99.9d0)
      r1 = min(r1,999.9d0)
      r2 = min(r2,999.9d0)
*
*       Print one line for every new stage (type change of #I1, #I2 or c.m.).
      WRITE (17,10)  tphys, name(i1), name(i2), kstar(i1), kstar(i2),
     &               kstar(i), m1, m2, r1, r2, ri, ecc, a0, p, iqcoll
   10 FORMAT (f8.1,2i6,3i4,2f5.1,f7.1,f6.1,f5.1,f7.3,f8.1,f9.1,i3)
      CALL flush(17)
*
   20 RETURN
*
      END