setup.f

Go to the documentation of this file.

      SUBROUTINE setup
*
*
*       Generation of initial coordinates & velocities.
*       -----------------------------------------------
*
      include 'common6.h'
      parameter(a0=1.0d0)
      REAL*8  ran2,a(8)
*
*
*       Choose between uniform density and Plummer model.
      kdum = idum1
      IF (kz(5).GT.0) go to 20
*
*       Set up a uniform spherical system.
      DO 10 i = 1,n
    1     a(1) = 0.0d0
          DO 2 k = 1,3
              a(k+1) = 2.0*ran2(kdum) - 1.0
              a(1) = a(1) + a(k+1)**2
    2     CONTINUE
          IF (a(1).GT.1.0) go to 1
*
    4     a(5) = 0.0d0
          DO 5 k = 1,3
              a(k+5) = 2.0*ran2(kdum) - 1.0
              a(5) = a(5) + a(k+5)**2
    5     CONTINUE
          IF (a(5).GT.1.0) go to 4
*
          DO 8 k = 1,3
*             X(K,I) = A(1)*A(K+1)
*       Density proportional to 1/R**2.
              x(k,i) = a(k+1)
*       Constant density.
              xdot(k,i) = a(k+5)
*       Isotropic velocities (magnitude randomized; no radial dependence).
    8     CONTINUE
   10 CONTINUE
*
      go to200 
*
*       Initialize centre of mass terms.
   20 DO 25 k = 1,3
          cmr(k) = 0.0d0
          cmrdot(k) = 0.0d0
   25 CONTINUE
*
*       Treat Plummer and lower options before Jaffe and Zhao models.
      IF (kz(5).GE.5) go to 52
*
*       Generate initial conditions from Plummer model (A & A 37, 183).
      DO 40 i = 1,n
   30     a(1) = ran2(kdum)
          IF (a(1).LT.1.0d-10) go to 30
          ri = (a(1)**(-0.6666667) - 1.0)**(-0.5)
*       Reject distant particles.
          IF (ri.GT.10.0) go to 30
*
          a(2) = ran2(kdum)
          a(3) = ran2(kdum)
          x(3,i) = (1.0 - 2.0*a(2))*ri
          x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
          x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
   32     a(4) = ran2(kdum)
          a(5) = ran2(kdum)
          a(6) = a(4)**2*(1.0 - a(4)**2)**3.5
          IF (0.1*a(5).GT.a(6)) go to 32
*
          a(8) = a(4)*sqrt(2.0)/(1.0 + ri**2)**0.25
          a(6) = ran2(kdum)
          a(7) = ran2(kdum)
          xdot(3,i) = (1.0 - 2.0*a(6))*a(8)
          xdot(1,i) = sqrt(a(8)**2 - xdot(3,i)**2)*cos(twopi*a(7))
          xdot(2,i) = sqrt(a(8)**2 - xdot(3,i)**2)*sin(twopi*a(7))
*
*       Accumulate c.m. terms.
          DO 35 k = 1,3
              cmr(k) = cmr(k) + body(i)*x(k,i)
              cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
   35     CONTINUE
   40 CONTINUE
*
*       Scale coordinates & velocities to analytical expectation values.
      sx = 1.5*twopi/16.0
      sv = sqrt(zmass/sx)
   42 DO 50 i = 1,n
          DO 45 k = 1,3
              x(k,i) = x(k,i) - cmr(k)/zmass
              xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
              x(k,i) = sx*x(k,i)
              xdot(k,i) = sv*xdot(k,i)
   45     CONTINUE
          IF (kz(22).EQ.1) THEN
              WRITE (10,48)  body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
   48         FORMAT (1p,7e14.6)
          END IF
   50 CONTINUE
      IF (kz(22).EQ.1) CALL flush(10)
*
*       Terminate standard cases (#5 <=1) or Jaffee/BH (#5 >= 5 label 52).
      IF (kz(5).LE.1.OR.kz(5).GE.5) go to 200
*
*       Save membership of first system for colour plot (N2 = NZERO - N1).
   52 n1 = n
*       Check optional generation of two orbiting Plummer models.
      IF (kz(5).EQ.2) THEN
          READ (5,*)  apo, ecc, n2, scale
          n2 = min(n,n2)
          semi = apo/(1.0 + ecc)
          semi = min(semi,50.0d0)
          semi = max(semi,2.0d0)
          ecc = min(ecc,0.999d0)
          ecc = max(ecc,0.0d0)
          zm2 = 0.0
          kskip = n1/n2
          DO 54 i = 1,n2
              j = kskip*i
              zm2 = zm2 + body(j)
   54     CONTINUE
          fac1 = zm2/(zmass + zm2)
          fac2 = zmass/(zmass + zm2)
*       Restrict volume ratio to 125 (i.e. unreasonable density contrast).
          IF (scale.LE.0.2d0) scale = 0.2
*       Increase total mass.
          zmass = zmass + zm2
*       Set apocentre velocity for new combined mass.
          vap = sqrt(zmass/semi)*sqrt((1.0 - ecc)/(1.0 + ecc))
          DO 55 i = 1,n
              IF (i.LE.n2) THEN
*       Copy members from first system by uniform skipping (N2 <= N1).
                  j = kskip*i
                  body(i+n) = body(j)
                  x(1,i+n) = scale*x(1,j) + fac2*apo
                  x(2,i+n) = scale*x(2,j)
                  x(3,i+n) = scale*x(3,j)
                  xdot(1,i+n) = xdot(1,j)/sqrt(scale)
                  xdot(2,i+n) = xdot(2,j)/sqrt(scale) + fac2*vap
                  xdot(3,i+n) = xdot(3,j)/sqrt(scale)
              END IF
              x(1,i) = x(1,i) - fac1*apo
              xdot(2,i) = xdot(2,i) - fac1*vap
   55     CONTINUE
      ELSE IF (kz(5).EQ.3) THEN
*       Prepare case of accretion disk with massive perturber.
          READ (5,*)  apo, ecc, dmin, scale
          rin = 0.5
          rout = 1.0
          zmass = 1.0
          body(1) = zmass
          DO 58 k = 1,3
              x(k,1) = 0.0
              xdot(k,1) = 0.0
   58     CONTINUE
*       Generate a thin disk population in circular orbits.
          DO 60 i = 2,n
              body(i) = 1.0d-03/float(n)
              semi = rin + (rout - rin)*float(i)/float(n)
              vcirc = sqrt((body(1) + body(i))/semi)
              phase = twopi*ran2(kdum)
              x(1,i) = semi*cos(phase)
              x(2,i) = semi*sin(phase)
              x(3,i) = 0.01*(2.0*ran2(kdum) - 1.0)
              xdot(1,i) = -vcirc*sin(phase)
              xdot(2,i) = vcirc*cos(phase)
              xdot(3,i) = 0.01*vcirc*(2.0*ran2(kdum) - 1.0)
   60     CONTINUE
*       Define membership of perturber and ensure no external tide.
          n2 = 1
          kz(14) = 0
*       Redefine solar mass unit and astronomical length scale in AU.
          zmbar = 1.0
          rbar = 1.0/2.05d+05
          body(n+1) = scale*body(1)
*       Set appropriate mass ratios for transforming to new c.m. frame.
          fac1 = body(n+1)/(zmass + body(n+1))
          fac2 = zmass/(zmass + body(n+1))
          zmass = zmass + body(n+1)
*       Form orbital elements for massive perturber (avoid ECC = 1).
          IF (abs(ecc - 1.0).GT.1.0d-05) THEN
              semi = dmin/(1.0 - ecc)
          ELSE
              semi = -1.0d+05
          END IF
          vm2 = zmass*(2.0/dmin - 1.0/semi)
          vap2 = zmass*(2.0/apo - 1.0/semi)
*       Determine initial y-velocity from angular momentum conservation.
          vy = sqrt(vm2)*dmin/apo
          vx = sqrt(vap2 - vy**2)
*       Place perturber on the Y-axis with appropriate velocities.
          x(1,n+1) = apo*fac2
          x(2,n+1) = 0.0
          x(3,n+1) = 0.0
          xdot(1,n+1) = -vx*fac2
          xdot(2,n+1) = vy*fac2
          xdot(3,n+1) = 0.0
*       Displace the disk members and include negative y-velocity.
          DO 70 i = 1,n
              x(1,i) = x(1,i) - fac1*apo
              xdot(1,i) = xdot(1,i) + fac1*vx
              xdot(2,i) = xdot(2,i) - fac1*vy
   70     CONTINUE
      ELSE IF (kz(5).EQ.4) THEN
*       Include two massive bodies (ECC > 1: NAME = 1 & 2 free floating).
          n2 = 0
          READ (5,*)  semi, ecc, zm1, zm2
          WRITE (6,72)  semi, ecc, zm1, zm2
   72     FORMAT (/,12x,'MASSIVE BODIES    A =',1p,e9.1,
     &            '  E =',0p,f6.2,'  M1/<M> =',f6.2,'  M2/<M> =',f6.2)
          body(1) = zm1
          body(2) = zm2
          IF (ecc.LT.1.0) THEN
*       Set apocentre velocity for new combined mass (using NAME = 1 & 2).
              vap = sqrt((zm1 + zm2)/semi)*sqrt((1.0 - ecc)/(1.0 + ecc))
              fac1 = zm2/(zm1 + zm2)
              fac2 = zm1/(zm1 + zm2)
              DO 75 k = 1,3
                  x(k,1) = 0.0
                  x(k,2) = 0.0
                  xdot(k,2) = 0.0
   75         CONTINUE
*       Initialize binary with c.m. at rest (elements change in SCALE).
              x(1,1) = -fac1*semi*(1.0 + ecc)
              x(1,2) = fac2*semi*(1.0 + ecc)
              xdot(2,1) = -fac1*vap
              xdot(2,2) = fac2*vap
          END IF
      ELSE IF (kz(5).EQ.5) THEN
*       Generate Jaffe model with length scale A0.
          zmh = 1.0/sqrt(float(2*n))
          DO 100 i = 1,n
   80         zm = ran2(kdum)
              ri = zm*a0/(1.0 - zm)
*       Reject distant particles.
              IF (ri.GT.10.0.OR.ri.LT.5.0d-03) go to 80
              a(2) = ran2(kdum)
              a(3) = ran2(kdum)
              x(3,i) = (1.0 - 2.0*a(2))*ri
              x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
              x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
              vi2 = 1.0/(ri + a0) + 2.0*zmh/ri
              DO 90 k = 1,3
                  vk2 = vi2*ran2(kdum)/3.0
                  zz = 2.0*ran2(kdum) - 1.0
                  ss = 1.0
                  IF (zz.LT.0.0) ss = -1.0
                  xdot(k,i) = ss*sqrt(vk2)
   90         CONTINUE
*       Accumulate c.m. terms.
              DO 95 k = 1,3
                  cmr(k) = cmr(k) + body(i)*x(k,i)
                  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
   95         CONTINUE
  100     CONTINUE 
          sx = 1.0
          sv = 1.0
*       Use existing procedure for c.m. correction.
          go to 42
      ELSE IF (kz(5).GE.6) THEN
*       Generate dwarf galaxy model with cusp & BH (Zhao, MN 278, 488).
          y0 = 5.0*1.25/twopi/sqrt(float(n))
*       Calculate RH and ZMH by setting mass inside RH to 5*N^{1/2}*<M>.
          rh = (exp(y0) - 1.0d0)**0.4
          zmh = 4.0*twopi/5.0*log(1.0d0 + rh**2.5)
*       Include possibility of using actual mass (read by SCALE).
          rcut = 0.0
          IF (kz(24).LT.0) THEN
              READ (5,*) zmh, rcut
              IF (zmh.GT.0.0) THEN
                  y0 = 1.25/twopi*(zmh/float(n))
                  rh = (exp(y0) - 1.0d0)**0.4
              END IF
          END IF
*
*       Add experimental choice of favourite value.
          IF (kz(28).EQ.3) THEN
              zmh = 1.0/sqrt(float(2*n))
              y0 = 1.25/twopi*(zmh/float(n))
              rh = (exp(y0) - 1.0d0)**0.4
          END IF
*
*       Define BH mass inside radius of influence (unscaled units).
          reject = 100.0
          cutm = 4.0*twopi/5.0*log(1.0d0 + reject**2.5)
*       Specify mass fraction using cut-off value and total in M_sun.
          zmh = zmh/(float(n)*zmbar)*cutm
          y0 = 1.25/twopi*zmh
          rh = (exp(y0) - 1.0d0)**0.4
          r2 = 100.0
          jdum = -1
          ranj = gasdev(jdum)
          IF (kz(5).EQ.7) THEN
              rh = 0.006
              zmh = log(1.0d0 + rh**0.75)
              reject = 200.0
              cutm = log(1.0d0 + reject**0.75)
              r2 = 200.0
          END IF
*
          icut = 0
          DO 120 i = 1,n
  110         zm = ran2(kdum)
              IF (kz(5).EQ.6) THEN
                  zx = 1.25/twopi*zm
                  zx = zx*cutm
                  zy = exp(zx)
                  ri = (zy - 1.0d0)**0.4
              ELSE
                  zx = zm*cutm
                  zy = exp(zx)
                  ri = (zy - 1.0d0)**(4.0/3.0)
              END IF
*       Reject large and small distances before scaling.
              IF (ri.GT.reject.OR.ri.LT.1.0d-04) go to 110
              a(2) = ran2(kdum)
              a(3) = ran2(kdum)
              x(3,i) = (1.0 - 2.0*a(2))*ri
              x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
              x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
  115         IF (kz(5).EQ.6) THEN
                  CALL vdisp(ri,r2,zmh,vi2)
              ELSE
                  CALL vdisp2(ri,r2,zmh,vi2)
              END IF
              sigma = sqrt(vi2)
              vx = sigma*gasdev(jdum)
              vy = sigma*gasdev(jdum)
              vz = sigma*gasdev(jdum)
              IF (max(abs(vx),abs(vy),abs(vz)).GT.3.0*sigma) go to 115
              xdot(1,i) = vx
              xdot(2,i) = vy
              xdot(3,i) = vz
              IF (ri.GT.rcut) go to 120
              semi = 2.0/ri - (vx**2 + vy**2 + vz**2)/zmh
              semi = 1.0/semi
              rrdot = x(1,i)*vx + x(2,i)*vy + x(3,i)*vz
              ecc2 = (1.0 - ri/semi)**2 + rrdot**2/(semi*zmh)
              ecc = sqrt(ecc2)
              IF (semi*(1.0 - ecc).LT.0.5*rcut) THEN
                  icut = icut + 1
                  go to 110
              END IF
  120     CONTINUE
*
          WRITE (6,125)  rh, zmh, cutm, rcut, icut
  125     FORMAT (/,12x,'BLACK HOLE    RH =',f6.3,'  MBH =',f8.4,
     &                '  CUTM =',f7.3,'  RCUT =',f9.5,'  ICUT =',i5)
*
          sx = 1.0
          sv = 1.0
*       Substitute c.m. value (hardly matters with zero velocity).
          IF (kz(24).LT.0) body(1) = zmh
          DO 130 k = 1,3
              x(k,1) = 0.0
              xdot(k,1) = 0.0
              cmr(k) = 0.0
              cmrdot(k) = 0.0
  130     CONTINUE
          DO 140 i = 1,n
              DO 135 k = 1,3
                  cmr(k) = cmr(k) + body(i)*x(k,i)
                  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
  135         CONTINUE
  140     CONTINUE
          go to 42
      END IF
*
*       Specify new membership.
      n = n + n2
      nzero = n
      ntot = n
      IF (n.GE.nmax-10) THEN
          WRITE (6,150)  n, nmax
  150     FORMAT (' DANGER!    LIMIT EXCEEDED   N =',i6,'  NMAX =',i6)
          stop
      END IF
*
      IF (kz(5).EQ.2) THEN
          WRITE (6,155)  semi, ecc, n1, n2, scale
  155     FORMAT (/,12x,'PLUMMER BINARY    A =',f6.2,'  E =',f6.2,
     &                  '  N1 =',i6,'  N2 =',i6,'  SCALE =',f6.2)
      ELSE IF (kz(5).EQ.3) THEN
          WRITE (6,160)  apo, ecc, dmin, scale
  160     FORMAT (/,12x,'MASSIVE PERTURBER    APO =',f6.2,'  E =',f6.2,
     &                  '  DMIN =',f6.2,'  MP/M1 =',f6.2)
      END IF
*
*       Re-initialize centre of mass terms.
      DO 162 k = 1,3
          cmr(k) = 0.0d0
          cmrdot(k) = 0.0d0
  162 CONTINUE
      zmass = 0.0
      DO 170 i = 1,n
          zmass = zmass + body(i)
          DO 165 k = 1,3
              cmr(k) = cmr(k) + body(i)*x(k,i)
              cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
  165     CONTINUE
  170 CONTINUE
      DO 180 i = 1,n
          DO 175 k = 1,3
              x(k,i) = x(k,i) - cmr(k)/zmass
              xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
  175     CONTINUE
  180 CONTINUE
*
*       Save random number sequence in COMMON for future use.
  200 idum1 = kdum
*
      RETURN
*
      END

      REAL*8 FUNCTION gasdev(IDUM)
C
C          GAUSSIAN DISTRIBUTION (PRESS P. 203).
C          -------------------------------------
C
      common/randg/  gset,iset
CC      DATA ISET/0/
      REAL*8 ran2
      IF (iset.EQ.0) THEN
1       v1=2.*ran2(idum)-1.
        v2=2.*ran2(idum)-1.
        r=v1**2+v2**2
        IF(r.GE.1.)go to 1
        fac=sqrt(-2.*log(r)/r)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      ELSE
        gasdev=gset
        iset=0
      ENDIF
      RETURN
      END

      SUBROUTINE vdisp(RI,R2,ZMH,VI2)
*
*       Velocity dispersion of dwarf galaxy model.
*       ------------------------------------------
*
      IMPLICIT REAL*8 (a-h,o-z)
*
*
      one3 = 1.0/3.0d0
      sum = 0.0
      rho = 1.0/(sqrt(ri)*(1.0 + ri**2*sqrt(ri)))
      ncall = 1000 + 100.0/ri**2
      dr = 2.0*(r2 - ri)/float(ncall)
      u = ri
*       Solve by Simpson's rule (NB! DR is twice the interval).
      DO 30 j = 1,ncall
      u1 = u
      CALL vfunc(zmh,u1,f1)
      u2 = u + 0.5d0*dr
      CALL vfunc(zmh,u2,f2)
      u3 = u + dr
      CALL vfunc(zmh,u3,f3)
      sum = sum + one3*(f1 + 4.0d0*f2 + f3)*dr/2.0d0
      u = u + dr
      IF (u.LT.1.0) THEN
          dr = 1.4*dr
      END IF
      IF (u.GT.100.0) go to 40
   30 CONTINUE
*
   40 vi2 = sum/rho
*
      RETURN
      END

      SUBROUTINE vfunc(ZMH,U,F)
      IMPLICIT REAL*8 (a-h,m,o-z)
*
      u2 = u**2
      zmu = 0.8*6.2831853*log(1.0d0 + u2*sqrt(u))
      rho = 1.0/(sqrt(u)*(1.0 + u2*sqrt(u)))
      f = rho*(zmu + zmh)/u2
*
      RETURN
      END

      SUBROUTINE vdisp2(RI,R2,ZMH,VI2)
*
*       Velocity dispersion of r^{-9/4} galaxy model.
*       ---------------------------------------------
*
      IMPLICIT REAL*8 (a-h,o-z)
*
*
      one3 = 1.0/3.0d0
      sum = 0.0
      rho = 1.0/(ri**2.25*(1.0 + ri**0.75))
      ncall = 1000 + 100.0/ri**2
      dr = 2.0*(r2 - ri)/float(ncall)
      u = ri
*       Solve by Simpson's rule (NB! DR is twice the interval).
      DO 30 j = 1,ncall
      u1 = u
      CALL vfunc2(zmh,u1,f1)
      u2 = u + 0.5d0*dr
      CALL vfunc2(zmh,u2,f2)
      u3 = u + dr
      CALL vfunc2(zmh,u3,f3)
      sum = sum + one3*(f1 + 4.0d0*f2 + f3)*dr/2.0d0
      u = u + dr
      IF (u.LT.1.0) THEN
          dr = 1.4*dr
      END IF
      IF (u.GT.200.0) go to 40
   30 CONTINUE
*
   40 vi2 = sum/rho
*
      RETURN
      END

      SUBROUTINE vfunc2(ZMH,U,F)
      IMPLICIT REAL*8 (a-h,m,o-z)
*
      u2 = u**2
      zmu = log(1.0d0 + u**0.75)
      rho = 1.0/(u**2.25*(1.0 + u**0.75))
      f = rho*(zmu + zmh)/u2
*
      RETURN
      END