setsys.f

Go to the documentation of this file.

      SUBROUTINE setsys
*
*
*       Selection of chain system.
*       --------------------------
*
      include 'common6.h'
      parameter(nmx=10,nmx3=3*nmx,nmxm=nmx*(nmx-1)/2)
      REAL*8  m,mass,mc,mij,mkk
      common/chain1/  xch(nmx3),vch(nmx3),m(nmx),
     &                zz(nmx3),wc(nmx3),mc(nmx),
     &                xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
     &                mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
      common/chainc/  xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
     &                listc(lmax)
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      common/binary/  zm(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)
      SAVE jsave
      INTEGER jsave(3)
*
*
*       Check whether new (or renewed) chain or addition of member(s).
      IF (nch.GT.0) go to 10
*
*       Include treatment for common envelope stage (denoted by NCH < 0).
      IF (nch.LT.0) THEN
          nch = -nch
*       Copy members from JLIST (set in CHTERM and modified in CMBODY).
          DO 1 l = 1,nch
              j = jlist(l)
              namec(l) = name(j)
              bodyc(l) = body(j)
              m(l) = body(j)
    1     CONTINUE
          go to 50
      END IF
*
*       Initialize chain indices, names & masses for maximum membership.
      DO 5 l = 1,4
          jlist(l) = 2*npairs + l
          namec(l) = name(2*npairs+l)
          bodyc(l) = body(2*npairs+l)
          m(l) = body(2*npairs+l)
    5 CONTINUE
      cm(9) = ebch0
*
*       Include treatment for near-synchronous binary as inert body (B-B).
      IF (jclose.GT.n.AND.kz(26).LT.2) THEN
          rsum = rmin
          nch = 2
          go to 10
      END IF
*
*       Define chain membership for three-body or four-body case.
      IF (jclose.LE.n.AND.jclose.GT.0) THEN
          nch = 3
          jlist(3) = jclose
          namec(3) = name(jclose)
          bodyc(3) = body(jclose)
          m(3) = body(jclose)
*       See whether a second single body should be added.
          IF (jcmax.LE.n.AND.jcmax.GE.ifirst) THEN
              jlist(4) = jcmax
              namec(4) = name(jcmax)
              bodyc(4) = body(jcmax)
              m(4) = body(jcmax)
              nch = 4
          END IF
      ELSE
          nch = 4
      END IF
*
*       Check for addition of binary (NCH <= 4).
      IF (jcmax.GT.n.AND.nch.LE.4) THEN
          ksp2 = jcmax - n
          IF (ksp2.GT.kspair) ksp2 = ksp2 - 1
          kspair = ksp2
          jcomp = 0
*       Save current members to prevent over-writing in KSTERM.
          DO 6 l = 1,nch
              jsave(l) = jlist(l)
    6     CONTINUE
          CALL ksterm
*       Note that second binary will now come first in N-body arrays.
          DO 7 l = 1,nch
              jlist(l) = jsave(l)
    7     CONTINUE
*       Add terminated KS components to chain arrays.
          DO 8 l = 1,2
              nch = nch + 1
              jlist(nch) = 2*npairs + l
              namec(nch) = name(2*npairs+l)
              bodyc(nch) = body(2*npairs+l)
              m(nch) = body(2*npairs+l)
    8     CONTINUE
      END IF
      go to 50
*
*       Improve coordinates & velocities of single perturber or c.m. body.
   10 CALL xvpred(jclose,-1)
*
*       Expand membership and save chain variables (single body or KS pair).
      IF (jclose.LE.n) THEN
          nch = nch + 1
          jlist(nch) = jclose
          namec(nch) = name(jclose)
          bodyc(nch) = body(jclose)
          m(nch) = body(jclose)
      ELSE
          kspair = jclose - n
*       Check for synchronous tidal binary (save dormant energy in ECOLL).
          IF (kz(27).GT.0.AND.kz(26).LT.2) THEN
              semi = -0.5*body(jclose)/h(kspair)
              IF (semi.LT.0.01*rsum) THEN
                  nch = nch + 1
                  jlist(nch) = jclose
                  namec(nch) = name(jclose)
                  bodyc(nch) = body(jclose)
                  m(nch) = body(jclose)
*       Define temporary KS ghost by saving masses and binding energy.
                  t0(2*kspair-1) = 1.0e+06
                  list(1,2*kspair-1) = 0
                  bodyc(9) = body(2*kspair-1)
                  bodyc(10) = body(2*kspair)
                  zmu = body(2*kspair-1)*body(2*kspair)/body(jclose)
                  ecoll = ecoll + zmu*h(kspair)
                  body(2*kspair-1) = 0.0d0
                  body(2*kspair) = 0.0d0
                  go to 50
              END IF
          END IF
*
*       Re-activate any merged binary before terminating as last pair.
          jg = 0
          IF (name(jclose).LT.0) THEN
*       Identify merger index and ghost for addition to chain.
              DO 12 k = 1,nmerge
                  IF (namem(k).EQ.name(jclose)) THEN
                      im = k
                  END IF
   12         CONTINUE
*       Note ghost must be single for maximum chain membership of 6.
              DO 14 j = 1,n
                  IF (body(j).EQ.0.0d0.AND.name(j).EQ.nameg(im)) THEN
                      jg = j
                  END IF
   14         CONTINUE
              WRITE (6,15)  name(jclose), name(jg), rsum,  r(jclose-n)
   15         FORMAT (' SETSYS HIARCH    NM NMG RSUM RB ',
     &                                   i6,i5,1p,2e10.2)
              iphase = 7
              CALL reset
              kspair = npairs
          END IF
*
*       Save global indices of existing members (KSTERM uses JLIST).
          DO 16 l = 1,nch
              jpert(l) = jlist(l)
   16     CONTINUE
*
*       Add energy of absorbed binary to the current initial energy.
          eb = body(2*kspair-1)*body(2*kspair)*h(kspair)/body(n+kspair)
          cm(9) = cm(9) + eb
          ebch0 = ebch0 + eb
*
*       Check saving of c.m. type and names.
          IF (ksave(1).EQ.0.OR.kstar(n+kspair).NE.0) THEN
              ksave(1) = kstar(n+kspair)
              ksave(2) = name(2*kspair-1) + name(2*kspair)
          END IF
*
*       Terminate KS pair and copy components (JCOMP=0 excludes ghost).
          iphase = 8
          jcomp = 0
          CALL ksterm
*
*       Restore old members.
          DO 18 l = 1,nch
              jlist(l) = jpert(l)
   18     CONTINUE
*
*       Add terminated KS components to chain arrays.
          DO 20 l = 1,2
              nch = nch + 1
              jlist(nch) = 2*npairs + l
              namec(nch) = name(2*npairs+l)
              bodyc(nch) = body(2*npairs+l)
              m(nch) = body(2*npairs+l)
   20     CONTINUE
*       See whether to include merger ghost.
          IF (jg.GT.0) THEN
              nch = nch + 1
              jlist(nch) = jg
              namec(nch) = name(jg)
              bodyc(nch) = body(jg)
              m(nch) = body(jg)
          END IF
          IF (nch.GT.6) THEN
              WRITE (6,30)  nch
   30         FORMAT (' DANGER!    NCH ',i4)
              stop
          END IF
      END IF
*
*       Specify membership for chain COMMON.
   50 nn = nch
*
      RETURN
*
      END