## 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
```

Stanford Exploration Project
11/18/1997