# SM implementation of algorithm used by Press et al. for deriving 
 # random numbers with gaussian distribution from random numbers with
 # uniform distribution, here the SM built-in function.  - bjw, December 2004

 # returns a vector of N randoms with rms=1.
 # gasdev <dimen> <random vector>
gasdev 2 #
  set dimen($2) = $1
  set _v1 = 2*random($1)-1
  set _v2 = 2*random($1)-1
  set _rsq = _v1**2 + _v2**2
  set _ifgood = (_rsq<1 && _rsq>0) ? 1 : 0
  set _fac = (_ifgood==1) ? (-2*ln(_rsq)/_rsq) : 0
  set _fac = sqrt(_fac)
  set _g1 = _v1 * _fac
  set _g2 = _v2 * _fac
  # try to fill up with good numbers from the two random vectors
  # if the entry wasn't good, we get a number from the second vector.
  # because more than half of the numbers fall in the unit circle,
  # this should work for any reasonably large N.
  define _j 0
  do _i = 0, $1-1 {
    if ( _ifgood[$_i]==1 ) { 
       set $2[$_i] = _g1[$_i] } \
    else { 
       while { _ifgood[$_j]==0 } {define _j ($_j+1) }
       set $2[$_i] = _g2[$_j] 
       define _j ($_j+1) }
  }


 # Copyright 2004, Benjamin Weiner, licensed under GPL.
 # You may use or redistribute this code freely, please retain this
 # notice.  If you incorporate it into a package, you can use the
 # package freely for your own research.  (You can't sell it without
 # distributing the source freely.)