      PROGRAM galaxy
*     model proposed by Schulman and Seiden
      IMPLICIT NONE
      REAL*8 p, s, two_pi, v
      REAL x, y
      INTEGER IR_, K, active_a, active_r, cell, dt, nactive, nring, t
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
      DATA K/0/
      CALL GWopen(IR_, 0)
      CALL parameter
      CALL initial
      DO WHILE(K .EQ. 0)
         CALL evolve
         CALL plot_spiral
         t = t + dt
         CALL GWevent(K, x, y)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE parameter
      IMPLICIT NONE
      REAL*8 pi, p, s, two_pi, v
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER active_a, active_r, cell, dt, nactive, nring, t
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
      nring = 50                            ! number of rings
      nactive = 200                         ! number of initial active cells
      v = 1                                 ! circular velocity
      p = 0.18D0                            ! star formation probability
      dt = 10                               ! time step
      s = 2*pi/6                            ! cell width
      two_pi = 2*pi
      END
C
      SUBROUTINE initial
      IMPLICIT NONE
      REAL*8 p, s, two_pi, v, x, y
      REAL xwin, ywin
      REAL dmy, rnd
      INTEGER active_a, active_r, cell, dt, nactive, nring,
     +       t, IR_, a, i, r, ICMB, Itmp1_, Itmp2_
      LOGICAL Ltmp1_, Ltmp2_
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
      COMMON /GWPLT/ICMB
      dmy = rnd(-1)
      DO Itmp1_ = 0, 50
         DO Itmp2_ = 1, 300
            cell(Itmp1_, Itmp2_) = 0
         END DO
      END DO
*     randomly activate cells
      i = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         Ltmp2_ = .TRUE.
         DO WHILE(Ltmp2_)
            x = int(nring*rnd(0)) + 1
            y = int(nring*rnd(0)) + 1
            r = int(sqrt(x*x + y*y)) + 1
            Ltmp2_ = r .GT. nring
         END DO
         a = int(6*r*rnd(0)) + 1            ! array index corresponds to an angle
         IF(cell(r,a) .EQ. 0) THEN
            i = i + 1
            active_a(i) = a                 ! location of active region
            active_r(i) = r
*           activate region, stars live for 15 time steps
            cell(r,a) = 15
         END IF
         Ltmp1_ = i .NE. nactive
      END DO
      t = 0                                 ! initial time
      CALL compute_aspect_ratio(REAL(nring), xwin, ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
*      CALL GWcolor(IR_, 19, 2)
*      CALL GWcolor(IR_, 0, 7)
      CALL GWcolor(IR_, 0, 2)
      CALL GWcolor(IR_, 19, 7)
      CALL GWsetpen(IR_, -8, -1, -1, 4)
      CALL GWsetbrs(IR_, -8, 1, -1)
      CALL GWbegincmb(IR_, 0, 1.0, 1.0)
      CALL GWellipse(IR_, 0.0,0.0,1.0,1.0)
      CALL GWendcmb(ICMB)
      END
C
      SUBROUTINE evolve
      IMPLICIT NONE
      REAL*8 angle, wt, p, s, two_pi, v
      INTEGER a, newactive_a(10000), newactive_r(10000), active_a,
     +        active_r, cell, dt, nactive, nring, t, am, ap, i,
     +        newactive, r, Itmp1_
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
*     number of active star clusters for next time step
      newactive = 0
      DO i = 1, nactive
         r = active_r(i)
         a = active_a(i)
*        activate neighboring cells of same ring
         CALL create(r,a+1,newactive,newactive_r,newactive_a)
         CALL create(r,a-1,newactive,newactive_r,newactive_a)
         angle = mod((a*s + v*t)/r, two_pi)
         IF(r .LT. nring) THEN
            wt = mod(v*t/(r+1), two_pi)
            ap = int((angle - wt)*(r+1)/s)
            ap = mod(ap,6*(r+1))
            CALL create(r+1,ap,newactive,newactive_r,newactive_a)
            CALL create(r+1,ap+1,newactive,newactive_r,newactive_a)
         END IF
         IF(r .GT. 1) THEN
            wt = mod(v*t/(r-1), two_pi)
            am = int((angle - wt)*(r-1)/s)
            am = mod(am, 6*(r-1))
            CALL create(r-1,am,newactive,newactive_r,newactive_a)
            CALL create(r-1,am+1,newactive,newactive_r,newactive_a)
         END IF
      END DO
      nactive = newactive
      DO Itmp1_ = 1, 10000
         active_r(Itmp1_) = newactive_r(Itmp1_)
      END DO
      DO Itmp1_ = 1, 10000
         active_a(Itmp1_) = newactive_a(Itmp1_)
      END DO
      END
C
      SUBROUTINE create(r,a,newactive,newactive_r,newactive_a)
      IMPLICIT NONE
      REAL*8 p, s, two_pi, v
      REAL rnd
      INTEGER newactive_r(10000), newactive_a(10000), active_a,
     +       active_r, cell, dt, nactive, nring, t, r, a, newactive
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
*     create star clusters
      IF(a .LT. 1) a = a + 6*r
      IF(rnd(0) .LT. p .AND. cell(r,a) .NE. 15) THEN
         newactive = newactive + 1
         newactive_a(newactive) = a
         newactive_r(newactive) = r
         cell(r,a) = 15                     ! activate cell
      END IF
      END
C
      SUBROUTINE plot_spiral
      IMPLICIT NONE
      CHARACTER Stmp1_*80
      REAL*8 p, s, two_pi, v
      REAL GWaspect, plotsize, theta, x, y
      INTEGER active_a, active_r, cell, dt,
     +       nactive, nring, t, IR_, a, r, ICMB
      COMMON /GLBLD/p, s, two_pi, v
      COMMON /GLBLI/active_a(10000), active_r(10000), cell(0:50,300),
     +             dt, nactive, nring, t
      COMMON /GWPLT/ICMB
      CALL GWsetogn(IR_, 0, -2)
      WRITE(Stmp1_,*) 'number of active star clusters = ', nactive
      CALL GWputtxt(IR_, -nring*GWaspect(1), REAL(nring), Stmp1_)
      DO r = 1, nring
         DO a = 1, 6*r
            IF(cell(r,a) .GT. 0) THEN
               theta = (a*s + v*t)/r
               x = r*cos(theta)
               y = r*sin(theta)
               plotsize = cell(r,a)/15.0
               CALL GWputcmb(IR_, ICMB, x, y, plotsize, plotsize, 0)
*              reduce star clusters lifetime
               cell(r,a) = cell(r,a) - 1
            END IF
         END DO
      END DO
      CALL GWflush(IR_, -2)
      END
