1      integer function ISRANP(XMEAN)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 2000-12-01 ISRANP Krogh  Removed unused parameter MONE.
6c>> 1996-03-30 ISRANP Krogh  Added external statement, rem. F90 comment.
7c>> 1995-11-16 ISRANP Krogh  Converted from SFTRAN to Fortran 77.
8c>> 1994-10-19 ISRANP Krogh  Changes to use M77CON
9c>> 1994-06-24 ISRANP CLL Changed common to use RANC[D/S]1 & RANC[D/S]2.
10c>> 1994-06-23 CLL Changed name to I[D/S]RANP.
11c>> 1992-03-16 CLL
12c>> 1991-11-26 CLL Reorganized common. Using RANCM[A/D/S].
13c>> 1991-11-22 CLL Added call to RAN0, and SGFLAG in common.
14c>> 1991-01-15 CLL Reordered common contents for efficiency.
15C>> 1990-01-23 CLL  Making names in common same in all subprogams.
16C>> 1987-04-23 IRANP  Lawson  Initial code.
17c        Returns one pseudorandom integer from the Poisson distribution
18c     with mean and variance = XMEAN.
19c        The probability of occurrence of the nonnegative integer k in
20c     the Poisson distribution with mean XMEAN is
21
22c                P(k) = exp(-XMEAN) * XMEAN**k / k!
23
24c        Let SUM(n) denote the sum for k = 0 to n of P(k).
25c     Let U be a random sample from a uniform
26c     distribution on [0.0, 1.0].  The returned value of ISRANP will be
27c     the smallest integer such that S(ISRANP) .ge. U.
28c     This variable has mean and variance = XMEAN.
29c     Reference: Richard H. Snow, Algorithm 342, Comm ACM, V. 11,
30c     Dec 1968, p. 819.
31c     Code based on subprogram written for JPL by Stephen L. Ritchie,
32c     Heliodyne Corp. and Wiley R. Bunton, JPL, 1969.
33c     Adapted to Fortran 77 for the JPL MATH77 library by C. L. Lawson &
34c     S. Y. Chiu, JPL, Apr 1987.
35c     ------------------------------------------------------------------
36c                    Subprogram argument and result
37c
38c     XMEAN [in]  Mean value for the desired distribution.
39c           Must be positive and not so large that exp(-XMEAN)
40c           will underflow.  For example, on a machine with underflow
41c           limit 0.14E-38, XMEAN must not exceed 88.
42c
43c     ISRANP [Returned function value]  Will be set to a nonnegative
44c           integer value if computation is successful.
45c           Will be set to -1 if XMEAN is out of range.
46c     ------------------------------------------------------------------
47c--S replaces "?": I?RANP, ?ERM1, ?ERV1, ?RANUA
48c--&                 RANC?1, RANC?2, ?PTR, ?NUMS, ?GFLAG
49c     Generic intrinsic functions referenced: EXP, LOG, MIN, NINT
50c     Other MATH77 subprogram referenced: RAN0
51c     Other MATH77 subprograms needed: ERMSG, ERFIN, AMACH
52c     Common referenced: /RANCS1/, /RANCS2/
53c     ------------------------------------------------------------------
54      external R1MACH
55      integer M, NMAX
56      parameter(M = 97, NMAX = 84)
57      real             ELIMIT, R1MACH
58      real             S, SNUMS(M), SPREV, SUM(0:NMAX)
59      real             TERM, TLAST, TWO
60      real             U, XMEAN, XMSAVE, ZERO
61      integer LAST, N, NMID
62      logical FIRST
63      parameter(ZERO = 0.0E0, TWO = 2.0E0)
64      common/RANCS2/SNUMS
65c
66      integer SPTR
67      logical SGFLAG
68      common/RANCS1/SPTR, SGFLAG
69c
70      save   /RANCS1/, /RANCS2/
71      save    ELIMIT, FIRST, LAST, NMID, SUM, TLAST, XMSAVE
72      data    FIRST / .true. /
73      data    ELIMIT / ZERO /
74      data    XMSAVE / 0.E0 /
75c     ------------------------------------------------------------------
76      if(FIRST .or. (XMEAN .ne. XMSAVE)) then
77c                                      Set for new XMEAN
78c     R1MACH(1) gives the underflow limit.
79c     ELIMIT gives a limit for arguments to the exponential function
80c     such that if XMEAN .le. ELIMIT,  exp(-XMEAN) should not
81c     underflow.
82c
83         if(ELIMIT .eq. ZERO) ELIMIT = -log(R1MACH(1) * TWO)
84         if(XMEAN .le. ZERO  .or.  XMEAN .gt. ELIMIT) then
85            call SERM1('ISRANP',1, 0,'Require 0 .lt. XMEAN .le. ELIMIT',
86     *        'XMEAN',XMEAN,',')
87            call SERV1('ELIMIT',ELIMIT,'.')
88            ISRANP = -1
89         else
90            LAST = 0
91            TLAST = exp(-XMEAN)
92            SUM(0) = TLAST
93            XMSAVE = XMEAN
94            NMID = nint(XMEAN)
95         end if
96         if(FIRST) then
97            FIRST = .false.
98            call RAN0
99         end if
100      end if
101c               SET U = UNIFORM RANDOM NUMBER
102      SPTR = SPTR - 1
103      if(SPTR .eq. 0) then
104         call SRANUA(SNUMS,M)
105         SPTR = M
106      endif
107      U = SNUMS(SPTR)
108c
109      if(U .gt. SUM(LAST)) then
110c                          COMPUTE MORE SUMS
111         N = LAST
112         SPREV = SUM(LAST)
113         TERM = TLAST
114   50    continue
115            N = N + 1
116            TERM = TERM * XMEAN / real(N)
117            S = SPREV + TERM
118            if(S .eq. SPREV) then
119               ISRANP = N
120               return
121            endif
122            if(N .le. NMAX) then
123               LAST = N
124               SUM(LAST) = S
125               TLAST = TERM
126            endif
127            if(U .le. S) then
128               ISRANP = N
129               return
130            endif
131            SPREV = S
132         go to 50
133      else
134c                            SEARCH THROUGH STORED SUMS
135c     Here we already know that U .le. SUM(LAST).
136c     It is most likely that U will be near SUM(NMID).
137c
138         if(NMID .lt. LAST) then
139            if( U .gt. SUM(NMID) ) then
140               do 100 N = NMID+1, LAST
141                  if(U .le. SUM(N)) then
142                     ISRANP = N
143                     return
144                  endif
145  100          continue
146            endif
147         endif
148c
149         do 150 N = min(NMID, LAST)-1, 0,-1
150            if(U .gt. SUM(N)) then
151               ISRANP = N+1
152               return
153            endif
154  150    continue
155         ISRANP = 0
156      endif
157      return
158      end
159