' ******************************************************************
' Gaussian random numbers
' ******************************************************************

DIM SHARED GaussSave AS DOUBLE  ' Saves a Gaussian number
DIM SHARED GaussNew AS INTEGER  ' Flags a new calculation

' ******************************************************************

FUNCTION RanGaussStd AS DOUBLE
' ------------------------------------------------------------------
' Computes 2 random numbers from the standard normal distribution,
' returns one and saves the other for the next call
' ------------------------------------------------------------------

  DIM AS DOUBLE X, R, Theta

  IF GaussNew = 0 THEN
    DO : X = RND : LOOP UNTIL X > 0
    R = SQR(-2 * LOG(X))

    DO : X = RND : LOOP UNTIL X > 0
    Theta = 6.283185307179586# * X

    RanGaussStd = R * COS(Theta)  ' Return 1st number
    GaussSave = R * SIN(Theta)    ' Save 2nd number
  ELSE
    RanGaussStd = GaussSave       ' Return saved number
  END IF

  GaussNew = NOT GaussNew
END FUNCTION

FUNCTION RanGauss(Mu AS DOUBLE, Sigma AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Returns a random number from a Gaussian distribution
' with mean Mu and standard deviation Sigma
' ------------------------------------------------------------------

  RanGauss = Mu + Sigma * RanGaussStd
END FUNCTION
