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