random_inv_gauss Function

FUNCTION random_inv_gauss(h, b, 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,INFINITY] FROM
  ! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
  ! WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
  ! USING A RATIO METHOD.
  !     H = PARAMETER OF DISTRIBUTION (0 <= REAL)
  !     B = PARAMETER OF DISTRIBUTION (0 < REAL)
  REAL, INTENT(IN)    :: h, b
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  !     Local variables
  REAL            :: ym, xm, r, w, r1, r2, x
  REAL, SAVE      :: a, c, d, e
  REAL, PARAMETER :: quart = 0.25
  IF (h < zero .OR. b <= zero) THEN
     WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     IF (h > quart*b*SQRT(vlarge)) THEN
        WRITE(*, *) 'THE RATIO H:B IS TOO SMALL'
        STOP
     END IF
     e = b*b
     d = h + one
     ym = (-d + SQRT(d*d + e))/b
     IF (ym < vsmall) THEN
        WRITE(*, *) 'THE VALUE OF B IS TOO SMALL'
        STOP
     END IF
     d = h - one
     xm = (d + SQRT(d*d + e))/b
     d = half*d
     e = -quart*b
     r = xm + one/xm
     w = xm*ym
     a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym))
     IF (a < vsmall) THEN
        WRITE(*, *) 'THE VALUE OF H IS TOO LARGE'
        STOP
     END IF
     c = -d*LOG(xm) - e*r
  END IF
  DO
     CALL RANDOM_NUMBER(r1)
     IF (r1 <= zero) CYCLE
     CALL RANDOM_NUMBER(r2)
     x = a*r2/r1
     IF (x <= zero) CYCLE
     IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT
  END DO
  fn_val = x
  RETURN
END FUNCTION random_inv_gauss