FUNCTION random_Poisson(mu, first) RESULT(ran_Poisson)
!**********************************************************************

!     Translated to Fortran 90 by Alan Miller from:
!                                    RANLIB
!
!           Library of Fortran Routines for Random Number Generation
!
!                           Compiled and Written by:
!
!                                Barry W. Brown
!                                 James Lovato
!
!                    Department of Biomathematics, Box 237
!                    The University of Texas, M.D. Anderson Cancer Center
!                    1515 Holcombe Boulevard
!                    Houston, TX      77030
!
! This work was supported by grant CA-16672 from the National Cancer Institute.

!                 GENerate POIsson random deviate Function

!  Generates a single random deviate from a Poisson
!  distribution with mean mu.

!                           Arguments

!  mu --> The mean of the Poisson distribution from which
!         a random deviate is to be generated.
!                           REAL mu

!                           Method

!  For details see:

!      Ahrens, J.H. and Dieter, U.
!      Computer Generation of Poisson Deviates
!      From Modified Normal Distributions.
!      ACM Trans. Math. Software, 8, 2 (June 1982), 163-179


!  MUPREV = PREVIOUS MU, MUOLD = MU AT LAST EXECUTION OF STEP P OR B.
!  TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
!  COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K) - DEL

!  SEPARATION OF CASES A AND B

IMPLICIT NONE

!     .. Scalar Arguments ..
REAL, INTENT(IN)    :: mu
LOGICAL, INTENT(IN) :: first
INTEGER             :: ran_Poisson

INTERFACE
  FUNCTION random_normal() RESULT (ran_norm)
    IMPLICIT NONE
    REAL :: ran_norm
  END FUNCTION random_normal
END INTERFACE

INTERFACE
  FUNCTION random_exponential() RESULT (ran_exp)
    IMPLICIT NONE
    REAL :: ran_exp
  END FUNCTION random_exponential
END INTERFACE

!     .. Local Scalars ..
REAL    :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,        &
           a4 = -.1661269, a5 = .1421878, a6 = -.1384794, a7 = .1250060,  &
           b1, b2, c, c0, c1, c2, c3, d, del, difmuk, e, fk, fx, fy, g,   &
           omega, p, p0, px, py, q, s, t, u, v, x, xx
INTEGER :: j, k, kflag, l, m
!     ..
!     .. Local Arrays ..
REAL    :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880. /)
REAL    :: pp(35)
!     ..
!     ..
!     .. Executable Statements ..
IF (.NOT. first) GO TO 10
IF (mu < 10.0) GO TO 120

!     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)

s = SQRT(mu)
d = 6.0*mu*mu

!     THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL PROBABILITIES FK
!     WHENEVER K >= M(MU).
!     L = INT(MU-1.1484) IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

l = INT(mu-1.1484)

!     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

10 g = mu + s*random_normal()
IF (g < 0.0) GO TO 20
ran_Poisson = INT(g)

!     STEP I. IMMEDIATE ACCEPTANCE IF ran_Poisson IS LARGE ENOUGH

IF (ran_Poisson.GE.l) RETURN

!     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U

fk = REAL(ran_Poisson)
difmuk = mu - fk
CALL RANDOM_NUMBER(u)
IF (d*u >= difmuk*difmuk*difmuk) RETURN

!     STEP P. PREPARATIONS FOR STEPS Q AND H.
!             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
!             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
!             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
!             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
!             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

20 IF (.NOT. first) GO TO 30
omega = .3989423/s
b1 = .4166667E-1 / mu
b2 = .3*b1*b1
c3 = .1428571*b1*b2
c2 = b2 - 15.*c3
c1 = b1 - 6.*b2 + 45.*c3
c0 = 1. - b1 + 3.*b2 - 15.*c3
c = .1069/mu
30 IF (g < 0.0) GO TO 50

!             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)

kflag = 0
GO TO 70

!     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN

!     STEP E. EXPONENTIAL SAMPLE - random_expon() FOR STANDARD EXPONENTIAL
!             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
!             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

50 e = random_exponential()
CALL RANDOM_NUMBER(u)
u = u + u - 1.0
t = 1.8 + SIGN(e, u)
IF (t <= (-.6744)) GO TO 50
ran_Poisson = INT(mu+s*t)
fk = REAL(ran_Poisson)
difmuk = mu - fk

!             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

kflag = 1
GO TO 70

!     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
RETURN

!     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
!             CASE ran_Poisson < 10 USES FACTORIALS FROM TABLE FACT

70 IF (ran_Poisson >= 10) GO TO 80
px = -mu
py = mu**ran_Poisson / fact(ran_Poisson + 1)
GO TO 110

!      CASE ran_Poisson >= 10 USES POLYNOMIAL APPROXIMATION
!      A0-A7 FOR ACCURACY WHEN ADVISABLE
!      .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

80 del = .8333333E-1 / fk
del = del - 4.8*del*del*del
v = difmuk/fk
IF (ABS(v) <= 0.25) GO TO 90
px = fk*LOG(1.0+v) - difmuk - del
GO TO 100

90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
100 py = .3989423/SQRT(fk)
110 x = (0.5-difmuk)/s
xx = x*x
fx = -0.5*xx
fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
IF (kflag > 0) THEN
  GO TO 60
ELSE
  GO TO 40
END IF

!     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)

120 IF (first) THEN
  m = MAX(1, INT(mu))
  l = 0
  p = EXP(-mu)
  q = p
  p0 = p
END IF

!     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

130 CALL RANDOM_NUMBER(u)
ran_Poisson = 0
IF (u <= p0) RETURN

!     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
!             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
!             (0.458 = PP(9) FOR MU = 10)

IF (l == 0) GO TO 150
j = 1
IF (u > 0.458) j = MIN(l, m)
DO k = j, l
  IF (u <= pp(k)) GO TO 180
END DO
IF (l == 35) GO TO 130

!     STEP C. CREATION OF NEW POISSON PROBABILITIES P
!             AND THEIR CUMULATIVES Q = PP(K)

150 l = l + 1
DO k = l, 35
  p = p*mu/REAL(k)
  q = q + p
  pp(k) = q
  IF (u <= q) GO TO 170
END DO
l = 35
GO TO 130

170 l = k
180 ran_Poisson = k
RETURN

END FUNCTION random_Poisson