program noise

! this program produces white noise with a selected rms value.
! the noise is very close to gaussian.

implicit real*8 (a - h, o - z)
character*40 output

write (*,'(5x, a)') 'enter filespec of output file.'
read (*, '(a)') output
open (unit = 1, file = output, status = 'unknown')

! input starting words for pseudo-random sequence (2 arbitrarily chosen
! integers); number of noise points to calculate; rms value of generated
! noise (it will approach this for a large number of generated points).

write (*,'(5x, a)') 'enter idum1, idum2, npoints, rms'
read (*, *) idum1, idum2, npoints, rms

idum1 = - abs (idum1)
idum2 = - abs (idum2)

sumnoise = 0.
do i = 1, npoints
  call  gauss0 (idum1, idum2, rms, xnoise)
  write (1, '(i6, 1pe13.5)') i, xnoise
  sumnoise = sumnoise + xnoise**2
end do
write (*, *) 'noise rms = ', sqrt (sumnoise / npoints)

close (unit = 1)
stop
end
!
subroutine gauss0 (idum1, idum2, sdev, x)

! generates normally distributed random number, zero-mean

! idum1 (input/output) must contain a negative number on the first entry
! idum2 (input/output) must contain a negative number on the first entry
! sdev (input) standard deviation of desired distribution
! x (output) the value of the computed normal variable

! box-muller method (polar algorithm):
! the algorithm is taken from donald e. knuth (standford university)
! "the art of computer programming", volume 2 / seminumerical algorithms,
! second edition, addison-wesley, isbn 0-201-03822-6, page 117

! 09-jun-88: first implementation

implicit real*8 (a - h, o - z)
save

10 v1 = 2.0 * (ran3 (idum1)) - 1.0
v2 = 2.0 * (ran3 (idum2)) - 1.0
s = v1 * v1 + v2 * v2
if (.not. (s .lt. 1.0)) go to 10
! select one of the two possible solutions
if (ran3 (idum1) .gt. 0.5) then
  x = sdev * v1 * dsqrt (- 2.0 * dlog (s) / s)
else
  x = sdev * v2 * dsqrt (- 2.0 * dlog (s) / s)
end if

return
end
!
function ran3 (idum)

implicit real*8 (a - h, o - z)
parameter (mbig = 1000000000, mseed = 161803398)
parameter (mz = 0, fac = 1. / mbig)
dimension ma (55)
data iff /0/
save

if (idum .lt. 0 .or. iff .eq. 0) then
  iff = 1
  mj = mseed - iabs (idum)
  mj = mod (mj, mbig)
  ma (55) = mj
  mk = 1
  do i = 1, 54
    ii = mod (21 * i, 55)
    ma (ii) = mk
    mk = mj - mk
    if (mk .lt. mz) mk = mk + mbig
    mj = ma (ii)
  end do
  do k = 1, 4
    do i = 1, 55
      ma (i) = ma (i) - ma (1 + mod (i + 30, 55))
      if (ma (i) .lt. mz) ma (i) = ma (i) + mbig
    end do
  end do
  inext = 0
  inextp = 31
  idum = 1
end if
inext = inext + 1
if (inext .eq. 56) inext = 1
inextp = inextp + 1
if (inextp .eq. 56) inextp = 1
mj = ma (inext) - ma (inextp)
if (mj .lt. mz) mj = mj + mbig
ma (inext) = mj
ran3 = mj * fac

return
end