Normal distribution

      SUBROUTINE sub5(npts,seed,output)

C        SUB5 generates gaussian deviates in OUTPUT.

C        independent
C        zero mean
C        unity variance

      INTEGER npts,seed
      REAL output(npts)

      INTEGER           j
      DOUBLE PRECISION  midpoint,scale,v1,v2,sumsq,squirt       
      DOUBLE PRECISION  dmagio1,dmagio2,dj1,dj2

      dmagio1 = DBLE(7*7*7*7*7)
      dmagio2 = DBLE(1 + 2*(2**30 - 1))
      dj2 = DBLE(seed)
      midpoint = DBLE(dmagio2)/2.0D0
      scale = 2.0D0/(DBLE(dmagio2) - 1.0D0)   
    
      DO j = 2,npts,2
111      dj1 = MOD(dmagio1*dj2,dmagio2)
         dj2 = MOD(dmagio1*dj1,dmagio2)

C           dj1 and dj2 cycle irregularly through the
C           integers 1 through (dmagio2 - 1)

         v1 = scale*(dj1 - midpoint)
         v2 = scale*(dj2 - midpoint)

C           v1 and v2 are the dj's scaled and shifted so that
C           they are uniformly distributed between -1 and 1

         sumsq = v1**2 + v2**2
         IF(sumsq.GE.1.0D0) THEN
            GO TO 111
         ELSE
         END IF

C           this IF tosses out all points {v1,v2} not inside
C           the unit circle

         squirt = SQRT((-2.0D0*LOG(sumsq))/sumsq)
         output(j - 1) = SNGL(v1*squirt)
         output(j) = SNGL(v2*squirt)
      END DO
      RETURN
      END