random_gamma1 Function

FUNCTION random_gamma1(s, first) RESULT(fn_val)
  ! Uses the algorithm in
  ! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
  ! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.
  ! Generates a random gamma deviate for shape parameter s >= 1.
  REAL, INTENT(IN)    :: s
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  ! Local variables
  REAL, SAVE  :: c, d
  REAL        :: u, v, x
  IF (first) THEN
     d = s - one/3.
     c = one/SQRT(9.0*d)
  END IF
  ! Start of main loop
  DO
     ! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.
     DO
        x = random_normal()
        v = (one + c*x)**3
        IF (v > zero) EXIT
     END DO
     ! Generate uniform variable U
     CALL RANDOM_NUMBER(u)
     IF (u < one - 0.0331*x**4) THEN
        fn_val = d*v
        EXIT
     ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN
        fn_val = d*v
        EXIT
     END IF
  END DO
  RETURN
END FUNCTION random_gamma1