hidat.f

Go to the documentation of this file.

      SUBROUTINE hidat
*
*
*       Hierarchical data bank.
*       -----------------------
*
      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),iflag(mmax)
      REAL*8  xx(3,3),vv(3,3),m1,m2,m3
*     REAL*4  EB(KMAX),ECM(KMAX)
      LOGICAL  first
      SAVE  first
      DATA  first /.true./
*
*
*       Write formatted data bank on unit 87.
      IF (first) THEN
          OPEN (unit=87,status='NEW',form='FORMATTED',file='HIDAT')
          first = .false.
          WRITE (87,1)
    1     FORMAT (/,'  NAM1  NAM2  NAM3  K*       M1   M2   M3    RI',
     &            '    EMAX   E0    E1      P0    P1')
      END IF
*
      IF (nmerge.GT.0) THEN
*       Count higher-order systems.
          mult = 0
          DO 2 i = n+1,ntot
              IF (name(i).LT.-2*nzero) mult = mult + 1
    2     CONTINUE
          WRITE (87,3)  npairs, nrun, n, nc, nmerge, mult, newhi, ttot
    3     FORMAT (/,i6,i4,i6,3i4,i6,f9.1)
      END IF
*
*       Produce output for each merged KS pair (TRIPLE, QUAD, .. [[B,B],B]).
      ipair = 0
*     ISTAB = 0
      DO 30 jpair = 1,npairs
          icm = n + jpair
          IF (name(icm).GT.0.AND.body(icm).GT.0.0d0) go to 30
          j2 = 2*jpair
          j1 = j2 - 1
*       Determine merger & ghost index (delay ghost c.m.).
          CALL findj(j1,j,im)
          kcm = kstarm(im)
          IF (body(icm).GT.0.0d0) THEN
              semi1 = -0.5*body(icm)/h(jpair)
              ecc1 = (1.0 - r(jpair)/semi1)**2 +
     &                                 tdot2(jpair)**2/(body(icm)*semi1)
          END IF
*
*       Consider any multiple hierarchies first (c.m. mass > 0 also here).
          IF (name(icm).LT.-2*nzero) THEN
              im = 0
*       Find the original merger and ghost index (J > N possible).
              DO 5 k = 1,nmerge
                  IF (namem(k).EQ.name(icm) + 2*nzero) im = k
    5         CONTINUE
              IF (im.EQ.0) go to 30
              j = 0
              DO 8 k = ifirst,ntot
                  IF (nameg(im).EQ.name(k)) j = k
    8         CONTINUE
              IF (j.EQ.0) go to 30
              ipair = ipair + 1
*       Employ actual masses and two-body distance for energy & eccentricity.
              bodycm = cm(1,im) + cm(2,im)
*             EB(IPAIR) = CM(1,IM)*CM(2,IM)*HM(IM)/BODYCM
              semi = -0.5*bodycm/hm(im)
              rj = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
              td2 = 0.0
              DO 10 k = 1,4
                  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
   10         CONTINUE
              ecc2 = (1.0 - rj/semi)**2 + td2**2/(bodycm*semi)
              m1 = cm(1,im)*smu
              m2 = cm(2,im)*smu
              m3 = body(j2)*smu
*       Determine EMAX and other relevant parameters.
              e0 = sqrt(ecc2)
              CALL himax(j1,im,e0,semi,emax,emin,zi,tg,edav)
              nam2 = name(j)
*       Check standard case of triple or quadruple (NAME(J2) <= N or > N).
          ELSE IF (body(icm).GT.0.0d0) THEN
              ipair = ipair + 1
*       Employ actual masses and two-body distance for energy & eccentricity.
              bodycm = cm(1,im) + cm(2,im)
*             EB(IPAIR) = CM(1,IM)*CM(2,IM)*HM(IM)/BODYCM
              semi = -0.5*bodycm/hm(im)
              rj = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
              td2 = 0.0
              DO 12 k = 1,4
                  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
   12         CONTINUE
              ecc2 = (1.0 - rj/semi)**2 + td2**2/(bodycm*semi)
*       Include separate diagnostics for the hierarchy (inner comps J1 & J).
              m1 = cm(1,im)*smu
              m2 = cm(2,im)*smu
              m3 = body(j2)*smu
*       Retain the next part just in case.
              e0 = sqrt(ecc2)
              e1 = sqrt(ecc1)
*       Determine EMAX and other relevant parameters (if needed).
              CALL himax(j1,im,e0,semi,emax,emin,zi,tg,edav)
              pmin = semi1*(1.0 - e1)
              IF (j.LT.0) j = j1
              rm = semi*(1.0 - e0)/max(radius(j1),radius(j),1.0d-20)
              rm = min(rm,99.9d0)
*       Obtain inclination between inner relative motion and outer orbit.
              DO 15 k = 1,3
                  xx(k,1) = xrel(k,im)
                  xx(k,2) = 0.0
                  xx(k,3) = x(k,j2)
                  vv(k,1) = vrel(k,im)
                  vv(k,2) = 0.0
                  vv(k,3) = xdot(k,j2)
   15         CONTINUE
              CALL inclin(xx,vv,x(1,icm),xdot(1,icm),angle)
              pcr = stability(m1,m2,m3,e0,e1,angle)*semi
              nam2 = name(j)
*       Perform stability check.
              IF (pmin*(1.0 - gamma(jpair)).LT.pcr.AND.e1.LT.0.96.AND.
     &            list(1,j1).GT.0) THEN
*                 ISTAB = JPAIR
                  WRITE (6,20)  name(j1), name(j), 180.0*angle/3.14, e1,
     &                          pmin, pcr, gamma(jpair), r0(jpair)
   20             FORMAT (' MERGE UNSTAB    NAM INC E1 PM PCR G R0 ',
     &                                      2i6,f7.1,f7.3,1p,4e10.2)
              END IF
          ELSE IF (body(j1).EQ.0.0d0) THEN
*       Treat case of ghost binary (JPAIR saved in CM(3,IM) & CM(4,IM).
              im = 0
*       Find the original merger index (use ghost name for [B,B]).
              DO 22 k = 1,nmerge
                  IF (nameg(k).EQ.name(icm)) im = k
   22         CONTINUE
              IF (im.EQ.0) go to 30
              kcm = kstar(icm)
              j = 0
*       Locate the first KS component (former c.m. hence subtract NZERO).
              DO 23 k = 1,ifirst
                  IF (name(k) - nzero.EQ.name(j1)) j = k
   23         CONTINUE
              IF (j.EQ.0) go to 30
*       Include the case of [[B,S],[B,S]] which requires more work.
              IF (body(j).EQ.0.0d0) THEN
                  jm = 0
                  DO 24 k = 1,nmerge
                      IF (namem(im).EQ.nameg(k)) jm = k
   24             CONTINUE
                  IF (jm.EQ.0) go to 30
                  j = 0
*       Employ new ghost identification to find the corresponding c.m index.
                  DO 25 k = n+1,ntot
                      IF (name(k).EQ.namem(jm)) j = k
   25             CONTINUE
                  IF (j.EQ.0) go to 30
                  j = 2*(j - n)
*       Note that both the triple masses (M1,M2) are saved in CM(1->4,JM).
              END IF
              ipair = ipair + 1
              nam2 = name(j2)
*       Set indices for the whole system (KS pair, c.m.; note NAME(J2) < 0).
              jp = kvec(j)
              icm = n + jp
              j2 = icm
              j = icm
*       Form outer semi-major axis and eccentricity.
              semi1 = -0.5*body(icm)/h(jp)
              ecc1 = (1.0 - r(jp)/semi1)**2 +
     &                                   tdot2(jp)**2/(body(icm)*semi1)
*       Copy masses from second merger component of [B,B] system.
              bodyj1 = cm(3,im)
              bodyj2 = cm(4,im)
              bodycm = bodyj1 + bodyj2
              bodycm = max(bodycm,1.0d-10)
              m1 = bodyj1*smu
              m2 = bodyj2*smu
              m3 = (body(icm) - bodycm)*smu
*             EB(IPAIR) = BODYJ1*BODYJ2*H(JPAIR)/BODYCM
              semi = -0.5*bodycm/h(jpair)
              ecc2 = (1.0 - r(jpair)/semi)**2 +
     &                                   tdot2(jpair)**2/(bodycm*semi)
              angle = 0.0
              emax = 0.0
          END IF
*
*       Evaluate the potential energy of c.m.
          phi = 0.0
          DO 26 k = ifirst,ntot
              IF (k.EQ.icm) go to 26
              rij2 = (x(1,k) - x(1,icm))**2 + (x(2,k) - x(2,icm))**2
     &                                      + (x(3,k) - x(3,icm))**2
              phi = phi + body(k)/sqrt(rij2)
   26     CONTINUE
*
*       Evaluate eccentricities and periods.
          e0 = sqrt(ecc2)
          e1 = sqrt(ecc1)
          p0 = days*semi*sqrt(abs(semi)/bodycm)
          p1 = days*semi1*sqrt(abs(semi1)/body(icm))
          p0 = min(p0,9999.0d0)
*       Obtain binding energy (per unit mass) of c.m. motion.
          vj2 = xdot(1,icm)**2 + xdot(2,icm)**2 + xdot(3,icm)**2
          IF (body(icm).EQ.0.0d0) vj2 = 0.0
*         ECM(IPAIR) = 0.5*VJ2 - PHI
          ri = sqrt((x(1,icm) - rdens(1))**2 + (x(2,icm) - rdens(2))**2
     &                                       + (x(3,icm) - rdens(3))**2)
          WRITE (87,28)  name(j1), nam2, name(j2), kstar(j1), kstar(j),
     &                   kcm, m1, m2, m3, ri, emax, e0, e1, p0, p1
   28     FORMAT (3i6,3i3,3f5.1,f7.2,f7.3,2f6.2,f8.1,1p,e9.1)
   30 CONTINUE
      CALL flush(87)
*
*       Check merger termination of unstable system.
*     IF (ISTAB.GT.0) THEN
*         KSPAIR = ISTAB
*         IPHASE = 7
*         CALL RESET
*     END IF
*
      RETURN
*
      END