cloud.f

Go to the documentation of this file.

      SUBROUTINE cloud(J)
*
*
*       Initialization of interstellar cloud.
*       -------------------------------------
*
      include 'common6.h'
      common/clouds/  xcl(3,mcl),xdotcl(3,mcl),bodycl(mcl),rcl2(mcl),
     &                clm(mcl),clmdot(mcl),cldot,vcl,sigma,rb2,pcl2,
     &                tcl,stepcl,ncl,newcl
      REAL*8  ran2,a(13)
*
*
*       Increase new cloud counter.
      newcl = newcl + 1
*       Generate a random position on the boundary sphere.
    1 a(1) = 0.0
      DO 2 k = 1,3
*       Note that IDUM1 is saved in COMMON6 for restarts.
          a(k+1) = 2.0*ran2(idum1) - 1.0
          a(1) = a(1) + a(k+1)**2
    2 CONTINUE
      IF (a(1).GT.1.0)  go to 1
      DO 3 k = 1,3
          xcl(k,j) = a(k+1)*sqrt(rb2/a(1))
    3 CONTINUE
      ranphi = twopi*ran2(idum1)
      randi = sqrt(ran2(idum1))
      randi = dasin(randi)
*       Azimuth & inclination for velocity from Henon (A & A, 19, 488).
      a(1) = sin(ranphi)
      a(2) = cos(ranphi)
      a(3) = sin(randi)
      a(4) = cos(randi)
      IF (kz(13).GE.2) THEN
          DO 10 k = 1,3
              a(5) = ran2(idum1)
              a(6) = twopi*ran2(idum1)
              a(k+10) = vcl + sigma*sqrt(-2.0d0*log(a(5)))*cos(a(6))
   10     CONTINUE
          a(10) = sqrt(a(11)**2 + a(12)**2 + a(13)**2)
      ELSE
          a(10) = vcl
      END IF
      a(11) = a(10)*a(2)*a(3)
      a(12) = a(10)*a(1)*a(3)
      a(13) = -a(10)*a(4)
      phi = datan2(-xcl(1,j),(-xcl(2,j)))
      a(1) = xcl(3,j)/sqrt(rb2)
      theta = dacos(a(1))
*       Eulerian angles of the transformation.
      a(1) = sin(phi)
      a(2) = cos(phi)
      a(3) = sin(theta)
      a(4) = cos(theta)
*       Generate isotropic velocities ignoring c.m. motion of cluster.
      xdotcl(1,j) = -a(4)*a(1)*a(11) + a(2)*a(12) - a(3)*a(1)*a(13)
      xdotcl(2,j) = -a(4)*a(2)*a(11) - a(1)*a(12) - a(3)*a(2)*a(13)
      xdotcl(3,j) = -a(3)*a(11) + a(4)*a(13)
*       Convert to rotating coordinates.
      xdotcl(1,j) = xdotcl(1,j) + 0.5*tidal(4)*xcl(2,j)
      xdotcl(2,j) = xdotcl(2,j) - 0.5*tidal(4)*xcl(1,j)
      bodycl(j) = 0.0
*       The cloud mass is increased smoothly from zero at boundary.
      DO 20 k = 1,3
          xcl(k,j) = xcl(k,j) - rdens(k)
   20 CONTINUE
*       Coordinates transformed from density centre to inertial frame.
*
      RETURN
*
      END