lagr.f

Go to the documentation of this file.

      SUBROUTINE lagr(C)
*
*
*       Lagrangian radii.
*       -----------------
*
      include 'common6.h'
      REAL*8 r2,mrt,rt
      common/work1/ r2(nmax)
      parameter(lx=11)
      REAL*8 c(3),flagr(lx),rlagr(lx),rm(lx),dens(lx),vr(lx),avm(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.9/
*
*
*       Set square radii of single particles & c.m. bodies.
      np = 0
      DO 10 i = ifirst,ntot
          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
   10 CONTINUE
*
*       Sort square distances with respect to the centre C.
      CALL sort1(np,r2,jlist)

*       Determine Lagrangian radii for specified mass fractions.
      zm = 0.0
      DO 20 il = 1,lx
          zm1 = zm
          zm = 0.0
          zmh = flagr(il)*zmass
          vm2 = 0.0
          IF (il.GT.1) THEN
              r1 = sqrt(r2(i))
              iprev = i
          ELSE
              r1 = 0.0
              iprev = 0
          END IF
          i = 0
   15     i = i + 1
          im = jlist(i)
          zm = zm + body(im)
*       Sum mass-weighted square velocity in current shell.
          IF (zm.GT.zm1) THEN
              vm2 = vm2 + body(im)*(xdot(1,im)**2 + xdot(2,im)**2 +
     &                                              xdot(3,im)**2)
          END IF
          IF (zm.LT.zmh) go to 15
          rlagr(il) = sqrt(r2(i))
          IF (kz(7).LT.5) go to 20
*       Form mean square velocity and density (print outer shell radii).
          dm = zm - zm1
          IF (abs(dm).LT.1.0d-10) THEN
              vr(il) = 0.0
          ELSE
              vr(il) = sqrt(vm2/dm)
          END IF
          dv = 2.0*twopi/3.0*(r2(i)**1.5 - r1**3)
          dens(il) = dm/dv
          IF (dens(il).LE.0.0d0) dens(il) = 1.0
          rm(il) = 0.5*(r1 + rlagr(il))
*       Obtain the average mass in solar units.
          IF (i.GT.iprev) THEN
              avm(il) = dm/float(i - iprev)*smu
          ELSE
              avm(il) = 0.0
          END IF
   20 CONTINUE
*
*       Obtain half-mass radius separately.
      zm = 0.0
      zmh = 0.5*zmass
      i = 0
   25 i = i + 1
      im = jlist(i)
      zm = zm + body(im)
      IF (zm.LT.zmh) go to 25
*
*       Replace approximate half-mass radius by actual value.
      rscale = sqrt(r2(i))
*
*       Determine tidal radius and mass if kz(14) = 2.
      if (kz(14).eq.2.and.tidal(1).ne.0.) then
      i = np
      mrt = zmass
   30 continue
      if (mrt.gt.0.) then
         rt = (mrt/tidal(1))**(1./3.)
      else
         rt = 0.0
      endif
      IF (sqrt(r2(i)).gt.rt.and.rt.gt.0.and.i.gt.1) then
         im = jlist(i)
         mrt = mrt - body(im)
         i = i - 1
         go to 30
      endif
      endif

      IF ((kz(7).EQ.2.OR.kz(7).EQ.4).AND.time.GE.tnext) THEN
          IF (kz(14).EQ.2) WRITE (6,*)  tphys,mrt,rt
      END IF
*
      IF (kz(7).GE.3.AND.time.GE.tnext) THEN
          IF (kz(14).EQ.2) WRITE (14,35)  tphys, mrt, rt
   35     FORMAT (3x,'TIDAL RADIUS    TPHYS MRT RT ',f8.1,1p,2e10.2)
      END IF
*
*       Check output options (line printer or unit 14 or both).
      IF (kz(7).EQ.2.OR.kz(7).EQ.4.AND.time.GE.tnext) THEN
          WRITE (6,40)  (log10(rlagr(k)),k=1,lx)
   40     FORMAT (/,' LAGR:  ',13f7.3)
      END IF
*
      IF (kz(7).GE.3.AND.time.GE.tnext) THEN
          WRITE (14,50)  ttot, (log10(rlagr(k)),k=1,lx)
   50     FORMAT ('  LAGR:  ',f7.1,13f7.3)
          CALL flush(14)
      END IF
*
      IF (kz(7).EQ.5.AND.time.GE.tnext) THEN
          WRITE (26,60)  ttot, (dens(k),k=1,lx)
   60     FORMAT ('  DENSITY (T =',f7.1,'): ',1p,13e10.2)
          WRITE (26,65)  (rm(k),k=1,lx)
   65     FORMAT ('  DISTANCE:  ',1p,13e10.2)
          CALL flush(26)
          WRITE (27,70)  ttot, (vr(k),k=1,lx)
   70     FORMAT ('  VELOCITY (T =',f7.1,'): ',1p,13e10.2)
          CALL flush(27)
          WRITE (36,75)  ttot, (avm(k),k=1,lx)
   75     FORMAT ('  AVERAGE MASS (T =',f7.1,'): ',13f7.3)
          CALL flush(36)
      END IF
      IF (kz(7).EQ.6.AND.time.GE.tnext) THEN
          WRITE (28,80)  ttot, (avm(k), rm(k),k=1,lx)
   80     FORMAT (' PAIRWISE <M> <RM> (T =',f7.1,'): ',
     &                                13(0p,f7.3,1p,e10.2))
          CALL flush(28)
      END IF
*
      RETURN
*
      END