core.f

Go to the documentation of this file.

      SUBROUTINE core
*
*
*       Density centre & core radius.
*       -----------------------------
*
      include 'common6.h'
      REAL*8  rlist(lmax),rho,rhos
      common/work1/  rho(nmax)
*
*
*       Set generous cutoff.
      rcore2 = rscale**2
      IF (time.GT.0.0d0) rcore2 = max(9.0d0*rc**2,rcore2)
*
*       Select > N/5 central particles.
    5 nc = 0
      DO 10 i = ifirst,ntot
          ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                   (x(3,i) - rdens(3))**2
          IF (ri2.LT.rcore2) THEN
              nc = nc + 1
              jlist(nc) = i
          END IF
   10 CONTINUE
*
      IF (nc.LT.n/5) THEN
          rcore2 = 1.5*rcore2
          go to 5
      END IF
*
*       Obtain individual densities.
      rhos = 0.0
      DO 50 l = 1,nc
          i = jlist(l)
          nnb = list(1,i)
*       Skip and set density to zero if there are no neighbours.
          IF (nnb.EQ.0) THEN
              rho(i) = 0.0
              go to 50
          END IF 
*       Estimate D6 to limit the candidates (additional factor of 2).
          d6 = 2.0*(6.0/float(nnb))**0.66*rs(i)**2
          xi = x(1,i)
          yi = x(2,i)
          zi = x(3,i)
   20     n6 = 0
*
          DO 25 lj = 1,nnb
              j = list(lj+1,i)
              rij2 = (xi - x(1,j))**2 + (yi - x(2,j))**2 +
     &                                  (zi - x(3,j))**2
              IF (rij2.LT.d6) THEN
                  n6 = n6 + 1
                  ilist(n6) = j
                  rlist(n6) = rij2
              END IF
   25     CONTINUE
*
          IF (n6.LT.6) THEN
*       Make another iteration in rare cases.
              d6 = 1.5*d6
              IF (n6.LT.nnb) go to 20
          END IF
*
*       Sort list of square distances.
          DO 40 ii = 1,n6
              DO 35 jj = ii+1,n6
                  IF (rlist(jj).LT.rlist(ii)) THEN
                      rdum = rlist(ii)
                      idum = ilist(ii)
                      rlist(ii) = rlist(jj)
                      ilist(ii) = ilist(jj)
                      rlist(jj) = rdum
                      ilist(jj) = idum
                  END IF
   35         CONTINUE
   40     CONTINUE
*
          i6 = min(n6,6)
*       Calculate total mass of five nearest neighbours.
          xmass = 0.0d0
          DO 45 lj = 1,i6-1
              xmass = xmass + body(ilist(lj))
   45     CONTINUE
*
          rho(i) = xmass/(rlist(i6)*sqrt(rlist(i6)))
*       Define particle number density (not used).
*         RHON = RHO(I)/ZMASS
*       Assign zero weight if not enough neighbours.
          IF (n6.LT.6) rho(i) = 0.0d0
          rhos = max(rhos,rho(i))
   50 CONTINUE
*
*       Determine the density centre.
      DO 60 k = 1,3
          rdens(k) = 0.0d0
   60 CONTINUE
*
      rho1 = 0.0d0
      DO 70 l = 1,nc
          i = jlist(l)
          DO 65 k = 1,3
              rdens(k) = rdens(k) + rho(i)*x(k,i)
   65     CONTINUE
          rho1 = rho1 + rho(i)
   70 CONTINUE
*
*       Set current density centre based on improved determination.
      rho1 = max(rho1,zmass/rscale**3)
      DO 75 k = 1,3
          rdens(k) = rdens(k)/rho1
*       Ignore density centre for Plummer binary models except at end.
          IF (kz(5).EQ.2.AND.ttot.LT.tcrit) THEN
              rdens(k) = 0.0
          END IF
   75 CONTINUE
*
*       Obtain density radius & averaged density.
      rc = 0.0d0
      rho2 = 0.0d0
      DO 80 l = 1,nc
          i = jlist(l)
          xid = x(1,i) - rdens(1)
          yid = x(2,i) - rdens(2)
          zid = x(3,i) - rdens(3)
          rid2 = xid**2 + yid**2 + zid**2
          rc = rc + rho(i)**2*rid2
          rho2 = rho2 + rho(i)**2
   80 CONTINUE
*
*       Form core radius and average & maximum density (scaled in OUTPUT).
      IF (rho2.GT.0.0d0) rc = sqrt(rc/rho2)
      rhod = rho2/rho1
      rhod = (3.0/12.566)*rhod
      rhom = (3.0/12.566)*rhos
*       NB! Scaling factor 5 of Casertano & Hut eq. (V.3) replaced by XMASS.
      IF (rc.EQ.0.0d0) THEN
          rc = rscale
          nc = n/2
      END IF
*
*       Set square core radius and its inverse for routine REGINT.
      rc2 = rc**2
      rc2in = 1.0/rc2
*
*       Sum particles & mass inside the core radius and set rms velocity.
      nc1 = 0
      vc = 0.d00
      zmc = 0.d00
      DO 85 l = 1,nc
          i = jlist(l)
          xid = x(1,i) - rdens(1)
          yid = x(2,i) - rdens(2)
          zid = x(3,i) - rdens(3)
          rid2 = xid**2 + yid**2 + zid**2
          IF (rid2.LT.rc2) THEN
              nc1 = nc1 + 1
              vc = vc + xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
              zmc = zmc + body(i)
          END IF
   85 CONTINUE
*
*       Set core membership & rms velocity.
      nc = max(nc1,2)
      vc = sqrt(vc/float(nc))
*
      RETURN
*
      END