random_beta Function

FUNCTION random_beta(aa, bb, first) RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
  ! FROM A BETA DISTRIBUTION WITH DENSITY
  ! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
  ! USING CHENG'S LOG LOGISTIC METHOD.
  !     AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
  !     BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
  REAL, INTENT(IN)    :: aa, bb
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  !     Local variables
  REAL, PARAMETER  :: aln4 = 1.3862944
  REAL             :: a, b, g, r, s, x, y, z
  REAL, SAVE       :: d, f, h, t, c
  LOGICAL, SAVE    :: swap
  IF (aa <= zero .OR. bb <= zero) THEN
     WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     a = aa
     b = bb
     swap = b > a
     IF (swap) THEN
        g = b
        b = a
        a = g
     END IF
     d = a/b
     f = a+b
     IF (b > one) THEN
        h = SQRT((two*a*b - f)/(f - two))
        t = one
     ELSE
        h = b
        t = one/(one + (a/(vlarge*b))**b)
     END IF
     c = a+h
  END IF
  DO
     CALL RANDOM_NUMBER(r)
     CALL RANDOM_NUMBER(x)
     s = r*r*x
     IF (r < vsmall .OR. s <= zero) CYCLE
     IF (r < t) THEN
        x = LOG(r/(one - r))/h
        y = d*EXP(x)
        z = c*x + f*LOG((one + d)/(one + y)) - aln4
        IF (s - one > z) THEN
           IF (s - s*z > one) CYCLE
           IF (LOG(s) > z) CYCLE
        END IF
        fn_val = y/(one + y)
     ELSE
        IF (4.0*s > (one + one/d)**f) CYCLE
        fn_val = one
     END IF
     EXIT
  END DO
  IF (swap) fn_val = one - fn_val
  RETURN
END FUNCTION random_beta