*
* $Id: poissn.F,v 1.1.1.1 1996/04/01 15:02:55 mclareni Exp $
*
* $Log: poissn.F,v $
* Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE POISSN (AMU,N,IERROR)
C
C    POISSON GENERATOR
C    CODED FROM LOS ALAMOS REPORT      LA-5061-MS
C    PROB(N)=EXP(-AMU)*AMU**N/FACT(N)
C        WHERE FACT(N) STANDS FOR FACTORIAL OF N
C    ON RETURN IERROR.EQ.0 NORMALLY
C              IERROR.EQ.1 IF AMU.LE.0.
C
      SAVE EXPMA,AMUOL,AMAX
      DATA AMUOL/-12345.67/
C    AMAX IS THE VALUE ABOVE WHICH THE NORMAL DISTRIBUTION MUST BE USED
#if defined(CERNLIB_IBM)
      DATA AMAX/170.0/
#endif
#if defined(CERNLIB_CRAY)
      DATA AMAX/5677./
#endif
#if defined(CERNLIB_VAX)||defined(CERNLIB_UNIX)
      DATA AMAX/88.0/
#endif
      IERROR= 0
      IF(AMU.GT.AMAX) GO TO 500
      IF(AMU.EQ.AMUOL) GO TO 200
      IF(AMU.GT.0.) GO TO 100
C    MEAN SHOULD BE POSITIVE
      IERROR=1
      N = 0
      GO TO 999
C    SAVE EXPONENTIAL FOR FURTHER IDENTICAL REQUESTS
  100 AMUOL=AMU
      EXPMA=EXP(-AMU)
  200 PIR=1.
      N=-1
  300 N=N+1
      PIR=PIR*RNDM(N)
      IF(PIR.GT.EXPMA) GO TO 300
      GO TO 999
C   NORMAL APPROXIMATION FOR AMU.GT.AMAX
  500 CALL RANNOR(RAN,B)
      N=RAN*SQRT(AMU)+AMU+.5
      GO TO 999
C   ENTRY FOR USER TO SET AMAX, SWITCHOVER POINT TO NORMAL APPROXIMATION
      ENTRY POISET(AMU)
      WRITE(6,1001) AMU
 1001 FORMAT(' POISSON RANDOM NUMBER GENERATOR TO SWITCH TO',
     +       ' NORMAL APPROXIMATION ABOVE AMU = ',F12.2)
      AMAX=AMU
 999  END