c	Function to return random integer from Poisson distn. with mean 'av'.
c	Uses random(is).	Is(2) is seed vector of length 2.
c
	function ipoiss (is, av)
	dimension is (2)
c
c
c
	dimension facl (0: 400)
	data facl(   0) /       0.000000000/
	data facl(   1) /       0.000000000/
	data facl(   2) /       0.693147181/
	data facl(   3) /       1.791759469/
	data facl(   4) /       3.178053830/
	data facl(   5) /       4.787491743/
	data facl(   6) /       6.579251212/
	data facl(   7) /       8.525161361/
	data facl(   8) /      10.604602903/
	data facl(   9) /      12.801827480/
	data facl(  10) /      15.104412573/
	data facl(  11) /      17.502307846/
	data facl(  12) /      19.987214496/
	data facl(  13) /      22.552163853/
	data facl(  14) /      25.191221183/
	data facl(  15) /      27.899271384/
	data facl(  16) /      30.671860106/
	data facl(  17) /      33.505073450/
	data facl(  18) /      36.395445208/
	data facl(  19) /      39.339884187/
	data facl(  20) /      42.335616461/
	data facl(  21) /      45.380138898/
	data facl(  22) /      48.471181352/
	data facl(  23) /      51.606675568/
	data facl(  24) /      54.784729398/
	data facl(  25) /      58.003605223/
	data facl(  26) /      61.261701761/
	data facl(  27) /      64.557538627/
	data facl(  28) /      67.889743137/
	data facl(  29) /      71.257038967/
	data facl(  30) /      74.658236349/
	data facl(  31) /      78.092223553/
	data facl(  32) /      81.557959456/
	data facl(  33) /      85.054467018/
	data facl(  34) /      88.580827542/
	data facl(  35) /      92.136175604/
	data facl(  36) /      95.719694542/
	data facl(  37) /      99.330612455/
	data facl(  38) /     102.968198615/
	data facl(  39) /     106.631760261/
	data facl(  40) /     110.320639715/
	data facl(  41) /     114.034211781/
	data facl(  42) /     117.771881400/
	data facl(  43) /     121.533081515/
	data facl(  44) /     125.317271149/
	data facl(  45) /     129.123933639/
	data facl(  46) /     132.952575036/
	data facl(  47) /     136.802722637/
	data facl(  48) /     140.673923648/
	data facl(  49) /     144.565743946/
	data facl(  50) /     148.477766952/
	data facl(  51) /     152.409592584/
	data facl(  52) /     156.360836303/
	data facl(  53) /     160.331128217/
	data facl(  54) /     164.320112263/
	data facl(  55) /     168.327445448/
	data facl(  56) /     172.352797139/
	data facl(  57) /     176.395848407/
	data facl(  58) /     180.456291418/
	data facl(  59) /     184.533828861/
	data facl(  60) /     188.628173424/
	data facl(  61) /     192.739047288/
	data facl(  62) /     196.866181673/
	data facl(  63) /     201.009316399/
	data facl(  64) /     205.168199483/
	data facl(  65) /     209.342586753/
	data facl(  66) /     213.532241495/
	data facl(  67) /     217.736934114/
	data facl(  68) /     221.956441819/
	data facl(  69) /     226.190548324/
	data facl(  70) /     230.439043566/
	data facl(  71) /     234.701723443/
	data facl(  72) /     238.978389562/
	data facl(  73) /     243.268849003/
	data facl(  74) /     247.572914096/
	data facl(  75) /     251.890402210/
	data facl(  76) /     256.221135550/
	data facl(  77) /     260.564940972/
	data facl(  78) /     264.921649799/
	data facl(  79) /     269.291097651/
	data facl(  80) /     273.673124286/
	data facl(  81) /     278.067573440/
	data facl(  82) /     282.474292688/
	data facl(  83) /     286.893133295/
	data facl(  84) /     291.323950094/
	data facl(  85) /     295.766601351/
	data facl(  86) /     300.220948647/
	data facl(  87) /     304.686856766/
	data facl(  88) /     309.164193580/
	data facl(  89) /     313.652829950/
	data facl(  90) /     318.152639620/
	data facl(  91) /     322.663499127/
	data facl(  92) /     327.185287704/
	data facl(  93) /     331.717887197/
	data facl(  94) /     336.261181979/
	data facl(  95) /     340.815058871/
	data facl(  96) /     345.379407062/
	data facl(  97) /     349.954118041/
	data facl(  98) /     354.539085519/
	data facl(  99) /     359.134205370/
	data facl( 100) /     363.739375556/
	data facl( 101) /     368.354496072/
	data facl( 102) /     372.979468886/
	data facl( 103) /     377.614197874/
	data facl( 104) /     382.258588773/
	data facl( 105) /     386.912549123/
	data facl( 106) /     391.575988217/
	data facl( 107) /     396.248817052/
	data facl( 108) /     400.930948279/
	data facl( 109) /     405.622296161/
	data facl( 110) /     410.322776527/
	data facl( 111) /     415.032306728/
	data facl( 112) /     419.750805600/
	data facl( 113) /     424.478193418/
	data facl( 114) /     429.214391867/
	data facl( 115) /     433.959323995/
	data facl( 116) /     438.712914186/
	data facl( 117) /     443.475088121/
	data facl( 118) /     448.245772745/
	data facl( 119) /     453.024896238/
	data facl( 120) /     457.812387981/
	data facl( 121) /     462.608178527/
	data facl( 122) /     467.412199572/
	data facl( 123) /     472.224383927/
	data facl( 124) /     477.044665493/
	data facl( 125) /     481.872979230/
	data facl( 126) /     486.709261137/
	data facl( 127) /     491.553448223/
	data facl( 128) /     496.405478487/
	data facl( 129) /     501.265290892/
	data facl( 130) /     506.132825342/
	data facl( 131) /     511.008022665/
	data facl( 132) /     515.890824588/
	data facl( 133) /     520.781173716/
	data facl( 134) /     525.679013516/
	data facl( 135) /     530.584288294/
	data facl( 136) /     535.496943180/
	data facl( 137) /     540.416924106/
	data facl( 138) /     545.344177791/
	data facl( 139) /     550.278651724/
	data facl( 140) /     555.220294147/
	data facl( 141) /     560.169054037/
	data facl( 142) /     565.124881095/
	data facl( 143) /     570.087725725/
	data facl( 144) /     575.057539025/
	data facl( 145) /     580.034272767/
	data facl( 146) /     585.017879389/
	data facl( 147) /     590.008311976/
	data facl( 148) /     595.005524249/
	data facl( 149) /     600.009470555/
	data facl( 150) /     605.020105849/
	data facl( 151) /     610.037385686/
	data facl( 152) /     615.061266207/
	data facl( 153) /     620.091704128/
	data facl( 154) /     625.128656731/
	data facl( 155) /     630.172081848/
	data facl( 156) /     635.221937855/
	data facl( 157) /     640.278183660/
	data facl( 158) /     645.340778693/
	data facl( 159) /     650.409682896/
	data facl( 160) /     655.484856711/
	data facl( 161) /     660.566261076/
	data facl( 162) /     665.653857411/
	data facl( 163) /     670.747607612/
	data facl( 164) /     675.847474040/
	data facl( 165) /     680.953419514/
	data facl( 166) /     686.065407302/
	data facl( 167) /     691.183401114/
	data facl( 168) /     696.307365094/
	data facl( 169) /     701.437263809/
	data facl( 170) /     706.573062246/
	data facl( 171) /     711.714725802/
	data facl( 172) /     716.862220279/
	data facl( 173) /     722.015511874/
	data facl( 174) /     727.174567173/
	data facl( 175) /     732.339353147/
	data facl( 176) /     737.509837142/
	data facl( 177) /     742.685986874/
	data facl( 178) /     747.867770425/
	data facl( 179) /     753.055156230/
	data facl( 180) /     758.248113081/
	data facl( 181) /     763.446610113/
	data facl( 182) /     768.650616800/
	data facl( 183) /     773.860102953/
	data facl( 184) /     779.075038710/
	data facl( 185) /     784.295394535/
	data facl( 186) /     789.521141209/
	data facl( 187) /     794.752249826/
	data facl( 188) /     799.988691789/
	data facl( 189) /     805.230438804/
	data facl( 190) /     810.477462876/
	data facl( 191) /     815.729736304/
	data facl( 192) /     820.987231676/
	data facl( 193) /     826.249921865/
	data facl( 194) /     831.517780024/
	data facl( 195) /     836.790779582/
	data facl( 196) /     842.068894242/
	data facl( 197) /     847.352097970/
	data facl( 198) /     852.640365001/
	data facl( 199) /     857.933669826/
	data facl( 200) /     863.231987192/
	data facl( 201) /     868.535292100/
	data facl( 202) /     873.843559798/
	data facl( 203) /     879.156765777/
	data facl( 204) /     884.474885771/
	data facl( 205) /     889.797895750/
	data facl( 206) /     895.125771919/
	data facl( 207) /     900.458490712/
	data facl( 208) /     905.796028792/
	data facl( 209) /     911.138363044/
	data facl( 210) /     916.485470574/
	data facl( 211) /     921.837328708/
	data facl( 212) /     927.193914982/
	data facl( 213) /     932.555207148/
	data facl( 214) /     937.921183163/
	data facl( 215) /     943.291821191/
	data facl( 216) /     948.667099599/
	data facl( 217) /     954.046996953/
	data facl( 218) /     959.431492015/
	data facl( 219) /     964.820563745/
	data facl( 220) /     970.214191292/
	data facl( 221) /     975.612353993/
	data facl( 222) /     981.015031375/
	data facl( 223) /     986.422203146/
	data facl( 224) /     991.833849198/
	data facl( 225) /     997.249949600/
	data facl( 226) /    1002.670484600/
	data facl( 227) /    1008.095434617/
	data facl( 228) /    1013.524780246/
	data facl( 229) /    1018.958502250/
	data facl( 230) /    1024.396581559/
	data facl( 231) /    1029.838999269/
	data facl( 232) /    1035.285736641/
	data facl( 233) /    1040.736775094/
	data facl( 234) /    1046.192096210/
	data facl( 235) /    1051.651681724/
	data facl( 236) /    1057.115513529/
	data facl( 237) /    1062.583573670/
	data facl( 238) /    1068.055844344/
	data facl( 239) /    1073.532307896/
	data facl( 240) /    1079.012946819/
	data facl( 241) /    1084.497743752/
	data facl( 242) /    1089.986681479/
	data facl( 243) /    1095.479742922/
	data facl( 244) /    1100.976911147/
	data facl( 245) /    1106.478169358/
	data facl( 246) /    1111.983500894/
	data facl( 247) /    1117.492889230/
	data facl( 248) /    1123.006317977/
	data facl( 249) /    1128.523770873/
	data facl( 250) /    1134.045231791/
	data facl( 251) /    1139.570684730/
	data facl( 252) /    1145.100113817/
	data facl( 253) /    1150.633503306/
	data facl( 254) /    1156.170837573/
	data facl( 255) /    1161.712101118/
	data facl( 256) /    1167.257278563/
	data facl( 257) /    1172.806354648/
	data facl( 258) /    1178.359314233/
	data facl( 259) /    1183.916142294/
	data facl( 260) /    1189.476823925/
	data facl( 261) /    1195.041344333/
	data facl( 262) /    1200.609688836/
	data facl( 263) /    1206.181842869/
	data facl( 264) /    1211.757791972/
	data facl( 265) /    1217.337521798/
	data facl( 266) /    1222.921018107/
	data facl( 267) /    1228.508266765/
	data facl( 268) /    1234.099253745/
	data facl( 269) /    1239.693965125/
	data facl( 270) /    1245.292387084/
	data facl( 271) /    1250.894505905/
	data facl( 272) /    1256.500307971/
	data facl( 273) /    1262.109779766/
	data facl( 274) /    1267.722907873/
	data facl( 275) /    1273.339678971/
	data facl( 276) /    1278.960079836/
	data facl( 277) /    1284.584097342/
	data facl( 278) /    1290.211718456/
	data facl( 279) /    1295.842930238/
	data facl( 280) /    1301.477719841/
	data facl( 281) /    1307.116074510/
	data facl( 282) /    1312.757981581/
	data facl( 283) /    1318.403428479/
	data facl( 284) /    1324.052402717/
	data facl( 285) /    1329.704891897/
	data facl( 286) /    1335.360883708/
	data facl( 287) /    1341.020365924/
	data facl( 288) /    1346.683326404/
	data facl( 289) /    1352.349753092/
	data facl( 290) /    1358.019634015/
	data facl( 291) /    1363.692957282/
	data facl( 292) /    1369.369711085/
	data facl( 293) /    1375.049883694/
	data facl( 294) /    1380.733463461/
	data facl( 295) /    1386.420438817/
	data facl( 296) /    1392.110798272/
	data facl( 297) /    1397.804530411/
	data facl( 298) /    1403.501623897/
	data facl( 299) /    1409.202067470/
	data facl( 300) /    1414.905849945/
	data facl( 301) /    1420.612960210/
	data facl( 302) /    1426.323387227/
	data facl( 303) /    1432.037120033/
	data facl( 304) /    1437.754147734/
	data facl( 305) /    1443.474459511/
	data facl( 306) /    1449.198044613/
	data facl( 307) /    1454.924892360/
	data facl( 308) /    1460.654992143/
	data facl( 309) /    1466.388333420/
	data facl( 310) /    1472.124905718/
	data facl( 311) /    1477.864698630/
	data facl( 312) /    1483.607701818/
	data facl( 313) /    1489.353905008/
	data facl( 314) /    1495.103297994/
	data facl( 315) /    1500.855870633/
	data facl( 316) /    1506.611612846/
	data facl( 317) /    1512.370514620/
	data facl( 318) /    1518.132566003/
	data facl( 319) /    1523.897757106/
	data facl( 320) /    1529.666078102/
	data facl( 321) /    1535.437519225/
	data facl( 322) /    1541.212070770/
	data facl( 323) /    1546.989723094/
	data facl( 324) /    1552.770466609/
	data facl( 325) /    1558.554291792/
	data facl( 326) /    1564.341189173/
	data facl( 327) /    1570.131149344/
	data facl( 328) /    1575.924162952/
	data facl( 329) /    1581.720220703/
	data facl( 330) /    1587.519313358/
	data facl( 331) /    1593.321431733/
	data facl( 332) /    1599.126566702/
	data facl( 333) /    1604.934709192/
	data facl( 334) /    1610.745850185/
	data facl( 335) /    1616.559980717/
	data facl( 336) /    1622.377091877/
	data facl( 337) /    1628.197174807/
	data facl( 338) /    1634.020220702/
	data facl( 339) /    1639.846220810/
	data facl( 340) /    1645.675166427/
	data facl( 341) /    1651.507048905/
	data facl( 342) /    1657.341859642/
	data facl( 343) /    1663.179590089/
	data facl( 344) /    1669.020231746/
	data facl( 345) /    1674.863776163/
	data facl( 346) /    1680.710214938/
	data facl( 347) /    1686.559539718/
	data facl( 348) /    1692.411742198/
	data facl( 349) /    1698.266814120/
	data facl( 350) /    1704.124747275/
	data facl( 351) /    1709.985533498/
	data facl( 352) /    1715.849164674/
	data facl( 353) /    1721.715632731/
	data facl( 354) /    1727.584929644/
	data facl( 355) /    1733.457047433/
	data facl( 356) /    1739.331978164/
	data facl( 357) /    1745.209713946/
	data facl( 358) /    1751.090246932/
	data facl( 359) /    1756.973569321/
	data facl( 360) /    1762.859673352/
	data facl( 361) /    1768.748551311/
	data facl( 362) /    1774.640195523/
	data facl( 363) /    1780.534598357/
	data facl( 364) /    1786.431752224/
	data facl( 365) /    1792.331649578/
	data facl( 366) /    1798.234282911/
	data facl( 367) /    1804.139644760/
	data facl( 368) /    1810.047727698/
	data facl( 369) /    1815.958524342/
	data facl( 370) /    1821.872027347/
	data facl( 371) /    1827.788229410/
	data facl( 372) /    1833.707123264/
	data facl( 373) /    1839.628701684/
	data facl( 374) /    1845.552957481/
	data facl( 375) /    1851.479883507/
	data facl( 376) /    1857.409472651/
	data facl( 377) /    1863.341717838/
	data facl( 378) /    1869.276612034/
	data facl( 379) /    1875.214148239/
	data facl( 380) /    1881.154319492/
	data facl( 381) /    1887.097118867/
	data facl( 382) /    1893.042539475/
	data facl( 383) /    1898.990574464/
	data facl( 384) /    1904.941217017/
	data facl( 385) /    1910.894460351/
	data facl( 386) /    1916.850297721/
	data facl( 387) /    1922.808722414/
	data facl( 388) /    1928.769727753/
	data facl( 389) /    1934.733307097/
	data facl( 390) /    1940.699453836/
	data facl( 391) /    1946.668161396/
	data facl( 392) /    1952.639423236/
	data facl( 393) /    1958.613232848/
	data facl( 394) /    1964.589583757/
	data facl( 395) /    1970.568469522/
	data facl( 396) /    1976.549883733/
	data facl( 397) /    1982.533820014/
	data facl( 398) /    1988.520272019/
	data facl( 399) /    1994.509233436/
	data facl( 400) /    2000.500697983/
c
c
c
c
	if (av .le. 0.0) then
		write (6,*) "From IPOISS: silly mean ",av
		ipoiss = -1
		return
		endif
c
c	Get int near peak of distn
	ip = av + 0.4999
c
c	If ip > 400, use a Normal approx
	if (ip .ge. 400) then
   10		continue
		ipoiss = anorm (is, av, sqrt(av)) + 0.4999
		if (ipoiss .lt. 0) goto 10
		return
		endif
c
c	Compute prob of getting value ip  as
c		exp(-av) * (av ** ip) / (ip !)
c	using table facl of log factorials
	pr = exp (-av + ip * alog(av) - facl(ip))
	pru = pr
	prd = pr
	ipu = ip
	ipd = ip
c
c	Get a uniform random
	z = random (is)
c
c	Accumulate probabilities of poisson values until z is exceeded
c	Values will be tried in order ip, ip-1, ip+1, ip-2, ip+2, etc.
c
	if (pr .lt. z) goto 20
	ipoiss = ip
	return
c
   20	continue
	if (ipd .eq. 0) goto 30
	prd = prd * ipd / av
	pr = pr + prd
	ipd = ipd - 1
	if (pr .lt. z) goto 30
	ipoiss = ipd
	return
c
   30	continue
	ipu = ipu + 1
	pru = pru * av / ipu
	pr = pr + pru
	if (pr .lt. z) goto 20
	ipoiss = ipu
	return
	end