lagr2.f

Go to the documentation of this file.

      SUBROUTINE lagr2(C)
*
*
*       Mass distribution for two mass groups.
*       --------------------------------------
*
      include 'common6.h'
      REAL*8 r2,zmi(nmax)
      common/work1/ r2(nmax)
      parameter(lx=11)
      REAL*8 c(3),flagr(lx),rlagr(lx)
*     DATA FLAGR/-1.9,-1.7,-1.5,-1.3,-1.1,-.9,-.7,-.5,-.3,-.1/
*     DATA FLAGR/0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,
*    &           0.75,0.9/
      DATA flagr/0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.625,0.75,0.99/
*
*
*       Set square radii of single particles & c.m. bodies (NAME <= NZERO/5).
      nam1 = nzero/5
*       Apply the mass test in same proportion for two Plummer spheres.
      IF (kz(5).EQ.2) THEN
          nam1 = (nzero - n1)/5
          nam2 = n1 + (nzero - n1)/5
      END IF
      iter = 0
      nblack = 0
      DO 5 i = 1,n
          IF (kstar(i).EQ.14) nblack = nblack + 1
    5 CONTINUE
*
    1 np = 0
      zm1 = 0.0
      DO 10 i = ifirst,ntot
          IF (kz(5).EQ.1) THEN
              IF (name(i).GT.nam1.AND.name(i).LE.nzero) go to 10
          ELSE IF (kz(5).EQ.2) THEN
              IF ((name(i).GT.nam1.AND.name(i).LE.n1).OR.
     &            (name(i).GT.nam2.AND.name(i).LE.nzero)) go to 10
          END IF
          zm1 = zm1 + body(i)
          np = np + 1
          r2(np) = (x(1,i) - c(1))**2 + (x(2,i) - c(2))**2 +
     &                                  (x(3,i) - c(3))**2
          jlist(np) = i
          ilist(np) = i
          zmi(np) = body(i)
   10 CONTINUE
*
*       Improve choice of NAM1 using relative deviation (max 10 tries).
      dm = (zm1 - 0.5*zmass)/zmass
      IF (abs(dm).GT.0.001.AND.iter.LE.10) THEN
          iter = iter + 1
          IF (abs(dm).GT.0.05) dm = dm*0.05/abs(dm)
          nm = dm*float(n-npairs)
          IF (nm.NE.0) THEN
              nam1 = nam1 - nm
              IF (kz(5).EQ.2) THEN
                  nam2 = nam2 - nm*(nzero - n1)/(5*nam2)
              END IF
              go to 1
          END IF
      END IF
*
*       Obtain sum of inverse separations to evaluate mass segregation.
      pot1 = 0.0
      DO 14 l = 1,np-1
          i = jlist(l)
          DO 12 ll = l+1,np
              j = jlist(ll)
              rij2 = 0.0
              DO 11 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
   11         CONTINUE
              pot1 = pot1 + 1.0/sqrt(rij2)
   12     CONTINUE
   14 CONTINUE
*       Define mean geometric distance.
      rp1 = float(np)*(np - 1)/(2.0*pot1)
*
*       Sort square distances with respect to the centre C.
      CALL sort1(np,r2,jlist)

*       Determine Lagrangian radii for specified mass fractions.
      DO 20 il = 1,lx
          zm = 0.0
          zmh = flagr(il)*zm1
          i = 0
   15     i = i + 1
          im = jlist(i)
          zm = zm + body(im)
          IF (zm.LT.zmh) go to 15
          rlagr(il) = sqrt(r2(i))
   20 CONTINUE
*
*       Obtain half-mass radius separately.
      zm = 0.0
      zmh = 0.5*zm1
      i = 0
   25 i = i + 1
      im = jlist(i)
      zm = zm + body(im)
      IF (zm.LT.zmh) go to 25

      rh1 = sqrt(r2(i))
      np1 = np
*
      WRITE (31,30)  time+toff, (log10(rlagr(k)),k=1,lx)
   30 FORMAT (' ',f8.1,13f7.3)
      CALL flush(31)
*
*       Sort masses from first group in increasing order.
      CALL sort1(np,zmi,ilist)
      zmx1 = zmi(1)
      zmx0 = 0.5*(zmx1 + body1)
      itry = 0
   32 np0 = 0
      lx0 = 0
*       Count dominant masses above half minimum + maximum single stars.
      DO 35 l = 1,np
          IF (zmi(l).GT.zmx0) THEN
              IF (lx0.EQ.0) lx0 = l
              np0 = np0 + 1
          END IF
   35 CONTINUE
*       Ensure at least 1% in the most massive bin or use NBLACK > 2.
      itry = itry + 1
      npcrit = sqrt(float(n - npairs))
      npcrit = min(npcrit,10)
      IF (nblack.GT.2) npcrit = nblack
      IF (np0.LT.npcrit.AND.itry.LT.20) THEN
          zmx0 = 0.9*zmx0
          go to 32
      END IF
      IF (np0.GT.npcrit.AND.itry.LT.20) THEN
          zmx0 = 1.1*zmx0
          go to 32
      END IF
*
*       Form sum of inverse separations to evaluate mass segregation.
      pot0 = 0.0
      lx0 = max(lx0,1)
      DO 44 l = lx0,np-1
          i = ilist(l)
          DO 42 ll = l+1,np
              j = ilist(ll)
              rij2 = 0.0
              DO 41 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
   41         CONTINUE
              pot0 = pot0 + 1.0/sqrt(rij2)
   42     CONTINUE
   44 CONTINUE
*       Define mean geometric distance.
      IF (np0.LE.1) pot0 = 1.0
      rp0 = float(np0)*(np0 - 1)/(2.0*pot0)
*
*       Search for single black holes and include KS binary.
      n14 = 0
      n14x = 0
      semi = 0
      potbh = 0.0
      zmbh = 0.0
      DO 45 i = ifirst,ntot
          IF (kstar(i).EQ.14.AND.body(i).GT.0.0d0) THEN
              n14 = n14 + 1
              ilist(n14) = i
              zmbh = zmbh + body(i)*smu
          ELSE IF (i.GT.n) THEN
*       Note: the following lines must be used for KS binary.
              j1 = 2*(i - n) - 1
              IF (kstar(j1).EQ.14.AND.kstar(j1+1).EQ.14) THEN
*       Add only one KS component to avoid small RIJ2 but use extra count.
                  n14 = n14 + 1
                  n14x = n14x + 1
                  ilist(n14) = i
                  zmbh = zmbh + body(i)*smu
                  ip = i - n
                  semi = -0.5*body(i)/h(ip)
              END IF
          END IF
   45 CONTINUE
*
*       Form sum of inverse separations to evaluate mass segregation.
      pot12 = 0.0
      DO 48 l = 1,n14-1
          i = ilist(l)
          DO 47 ll = l+1,n14
              j = ilist(ll)
              rij2 = 0.0
              DO 46 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
   46         CONTINUE
              potbh = potbh + 1.0/sqrt(rij2)
*       Note internal terms not included at present.
              IF (i.LT.ifirst) pot12 = pot12 + 1.0/sqrt(rij2)
   47     CONTINUE
   48 CONTINUE
*       Define mean geometric distance (modified by binary 10/2/11).
      n14 = n14 + n14x
      IF (n14.GT.1) rbh = float(n14)*(n14 - 1)/(2.0*potbh)
*       Omit internal binary contribution for secondary definition.
      IF (pot12.GT.0.0) THEN
          rx = float(n14-1)*(n14 - 2)/(2.0*(potbh-pot12))
      ELSE
          rx = 0.0
      END IF
*
*       Search for single neutron stars and include KS binary.
      n13 = 0
      semi2 = 0
      pot13 = 0.0
      r13 = 0.0
      zm13 = 0.0
      DO 145 i = ifirst,ntot
          IF (kstar(i).EQ.13) THEN
              n13 = n13 + 1
              ilist(n13) = i
              zm13 = zm13 + body(i)*smu
          ELSE IF (i.GT.n) THEN
              j1 = 2*(i - n) - 1
              IF (kstar(j1).EQ.13.AND.kstar(j1+1).EQ.13) THEN
                  n13 = n13 + 1
                  ilist(n13) = i
                  zm13 = zm13 + body(i)*smu
                  ip = i - n
                  semi2 = -0.5*body(i)/h(ip)
              END IF
          END IF
  145 CONTINUE
*
*       Form sum of inverse separations to evaluate mass segregation.
      DO 148 l = 1,n13-1
          i = ilist(l)
          DO 147 ll = l+1,n13
              j = ilist(ll)
              rij2 = 0.0
              DO 146 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
  146         CONTINUE
              pot13 = pot13 + 1.0/sqrt(rij2)
  147     CONTINUE
  148 CONTINUE
*       Define mean geometric distance (modified by binary 10/2/11).
      IF (n13.GT.1) r13 = float(n13)*(n13 - 1)/(2.0*pot13)
*
*       Treat the low-mass particles in the same way (exclude c.m. bodies).
      np = 0
      zm2 = 0.0
      DO 50 i = ifirst,n
          IF (kz(5).EQ.1) THEN
              IF (name(i).LE.nam1) go to 50
          ELSE IF (kz(5).EQ.2) THEN
              IF (name(i).LE.nam1.OR.
     &           (name(i).GT.n1.AND.name(i).LE.nam2)) go to 50
          END IF
          zm2 = zm2 + body(i)
          np = np + 1
          r2(np) = (x(1,i) - c(1))**2 + (x(2,i) - c(2))**2 +
     &                                  (x(3,i) - c(3))**2
          jlist(np) = i
   50 CONTINUE
*
*       Treat second mass group in the same way.
      pot2 = 0.0
      DO 54 l = 1,np-1
          i = jlist(l)
          DO 52 ll = l+1,np
              j = jlist(ll)
              rij2 = 0.0
              DO 51 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
   51         CONTINUE
              pot2 = pot2 + 1.0/sqrt(rij2)
   52     CONTINUE
   54 CONTINUE
      rp2 = float(np)*(np - 1)/(2.0*pot2)
*
*       Sort square distances with respect to the centre C.
      CALL sort1(np,r2,jlist)

      DO 60 il = 1,lx
          zm = 0.0
          zmh = flagr(il)*zm2
          i = 0
   55     i = i + 1
          im = jlist(i)
          zm = zm + body(im)
          IF (zm.LT.zmh) go to 55
          rlagr(il) = sqrt(r2(i))
   60 CONTINUE
*
*       Obtain half-mass radius separately.
      zm = 0.0
      zmh = 0.5*zm2
      i = 0
   65 i = i + 1
      im = jlist(i)
      zm = zm + body(im)
      IF (zm.LT.zmh) go to 65
*
      rh2 = sqrt(r2(i))
      np2 = np
*
      WRITE (32,30)  time+toff, (log10(rlagr(k)),k=1,lx)
      CALL flush(32)
*
      WRITE (6,70)  time+toff, np0, np1, np2, rh1, rh2, rp0, rp1, rp2
   70 FORMAT(/,' MASS GROUPS:    T NP0 NP1 NP2 RM1 RM2 RP0 RP1 RP2 ',
     &                           f7.1,i5,i6,i7,1x,5f7.3)
*
      IF (n14.GT.1.OR.n13.GT.1) THEN
          n14 = max(n14,1)
          WRITE (6,80)  time+toff, tphys, n14, zmbh/float(n14), rbh,
     &                  semi, n13, r13, rx
   80     FORMAT (/,' BH/NS SUBSYSTEM:    T TPHYS NBH <MBH> RBH A ',
     &                                    'NS R13 RX ',
     &                       2f7.1,i4,f6.1,f7.3,1p,e10.2,0p,i5,2f7.3)
      END IF
*
      RETURN
*
      END