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