]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:55 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE POISSN (AMU,N,IERROR) | |
11 | C | |
12 | C POISSON GENERATOR | |
13 | C CODED FROM LOS ALAMOS REPORT LA-5061-MS | |
14 | C PROB(N)=EXP(-AMU)*AMU**N/FACT(N) | |
15 | C WHERE FACT(N) STANDS FOR FACTORIAL OF N | |
16 | C ON RETURN IERROR.EQ.0 NORMALLY | |
17 | C IERROR.EQ.1 IF AMU.LE.0. | |
18 | C | |
19 | SAVE EXPMA,AMUOL,AMAX | |
20 | DATA AMUOL/-12345.67/ | |
21 | C AMAX IS THE VALUE ABOVE WHICH THE NORMAL DISTRIBUTION MUST BE USED | |
22 | #if defined(CERNLIB_IBM) | |
23 | DATA AMAX/170.0/ | |
24 | #endif | |
25 | #if defined(CERNLIB_CRAY) | |
26 | DATA AMAX/5677./ | |
27 | #endif | |
28 | #if defined(CERNLIB_VAX)||defined(CERNLIB_UNIX) | |
29 | DATA AMAX/88.0/ | |
30 | #endif | |
31 | IERROR= 0 | |
32 | IF(AMU.GT.AMAX) GO TO 500 | |
33 | IF(AMU.EQ.AMUOL) GO TO 200 | |
34 | IF(AMU.GT.0.) GO TO 100 | |
35 | C MEAN SHOULD BE POSITIVE | |
36 | IERROR=1 | |
37 | N = 0 | |
38 | GO TO 999 | |
39 | C SAVE EXPONENTIAL FOR FURTHER IDENTICAL REQUESTS | |
40 | 100 AMUOL=AMU | |
41 | EXPMA=EXP(-AMU) | |
42 | 200 PIR=1. | |
43 | N=-1 | |
44 | 300 N=N+1 | |
45 | PIR=PIR*RNDM(N) | |
46 | IF(PIR.GT.EXPMA) GO TO 300 | |
47 | GO TO 999 | |
48 | C NORMAL APPROXIMATION FOR AMU.GT.AMAX | |
49 | 500 CALL RANNOR(RAN,B) | |
50 | N=RAN*SQRT(AMU)+AMU+.5 | |
51 | GO TO 999 | |
52 | C ENTRY FOR USER TO SET AMAX, SWITCHOVER POINT TO NORMAL APPROXIMATION | |
53 | ENTRY POISET(AMU) | |
54 | WRITE(6,1001) AMU | |
55 | 1001 FORMAT(' POISSON RANDOM NUMBER GENERATOR TO SWITCH TO', | |
56 | + ' NORMAL APPROXIMATION ABOVE AMU = ',F12.2) | |
57 | AMAX=AMU | |
58 | 999 END |