PROGRAM TESTY
C
C  Test the poisson random number generator
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 CNT
      PARAMETER ( CNT = 100 )
      REAL *4 DATA(0:CNT), PROB(0:CNT), XPLOT(0:CNT), DELTA(0:CNT)
      REAL *4 XMIN, XMAX, YMIN, YMAX, XHI, XLO, YHI, YLO, MEAN
      CHARACTER *40 TITLE
      REAL *8 LAMBDA, VALUE
      INTEGER *4 GETCOUNT, NPTS, I, IX
      EXTERNAL GETCOUNT
  100 CONTINUE
      WRITE(6,*) ' Input mean and number of samples '
      READ (5,*) MEAN, NPTS
      IF ( NPTS .LE. 0 ) GO TO 900
      LAMBDA = MEAN
      DO I = 0, CNT
          DATA(I) = 0.
          CALL POISSON ( I, LAMBDA, VALUE )
          PROB(I) = FLOAT(NPTS) * VALUE
          DELTA(I) = -PROB(I)
          XPLOT(I) = FLOAT(I)
      END DO
      DO I = 1, NPTS
          IX = GETCOUNT ( MEAN )
          IX = MIN ( CNT, MAX( 0, IX ) )
          DATA(IX) = DATA(IX) + 1
          DELTA(IX) = 10.*( DATA(IX) - PROB(IX) )
      END DO
      YMIN = 0.0
      YMAX = 0.0
      XMIN = CNT
      XMAX = 0.0
      DO I = 1, CNT
          YMAX = MAX ( YMAX, DATA(I) )
          YMAX = MAX ( YMAX, PROB(I) )
          YMAX = MAX ( YMAX, DELTA(I) )
          YMIN = MIN ( YMIN, DELTA(I) )
          IF ( PROB(I) .GT. 0.0 ) THEN
              XMAX = MAX ( XMAX, FLOAT(I) )
              XMIN = MIN ( XMIN, FLOAT(I) )
          END IF
      END DO
      DO I = 0, CNT
          WRITE(7,1200) I, PROB(I), DATA(I), DATA(I)-PROB(I)
      END DO
C
C    Make the plot
C
      CALL PGBEGIN (0, '?', 1, 1 )
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
      CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
      CALL PGENV   ( XLO, XHI, YLO, YHI, 0, 2 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'POISSON DISTRIBUTION')
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, '(X)')
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, '(Y)')
      WRITE(TITLE,'(I5,A)')  NPTS, ' SAMPLES '
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE )
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, ' ' )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, ' ' )
      CALL PGSCI   ( 2 )
      CALL PGLINE  ( NPTS, XPLOT(0), PROB(0) )
      CALL PGPOINT ( NPTS, XPLOT(0), DATA(0), 4 )
      CALL PGEND
      GO TO 100
  900 CONTINUE
      STOP
 1200 FORMAT ( I5, 3F10.3 )
      END
      INTEGER *4 FUNCTION GETCOUNT ( MEAN )
C
C   Returns a poisson distributed random integer of mean MEAN
C   The poisson distribution is truncated at SIGMA standard deviations
C   above the mean.
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *4 MEAN, SIGMA, RAN0, LOLIMIT, HILIMIT, Y
      REAL *8 LAMBDA, PROB
      INTEGER *4 IX, ISEED
      EXTERNAL RAN0
      DATA ISEED, SIGMA  / 34297, 5.  /
      IF ( MEAN .LE. 0. ) THEN
          GETCOUNT = 0
          GO TO 900
      END IF
      LOLIMIT  = 0.
      HILIMIT  = MEAN + SIGMA*SQRT(MEAN)
      LAMBDA = MEAN
  100 CONTINUE
          IX = INT ( LOLIMIT + (HILIMIT-LOLIMIT) * RAN0(ISEED) )
          Y  = RAN0(ISEED)
          CALL POISSON ( IX, LAMBDA, PROB )
      IF ( PROB .LT. Y ) GO TO 100
      GETCOUNT = IX
  900 CONTINUE
      RETURN
      END
      SUBROUTINE POISSON ( IX, LAMBDA, PROB )
C
C   PROB is the probability for a random integer variable,
C   IX,   distributed in a poisson distribution of mean LAMBDA.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 IX, I
      REAL *8 LAMBDA, PROB, X, PI, FACT(0:50)
      LOGICAL FIRST
      SAVE
      DATA FIRST / .TRUE. /
      IF ( FIRST ) THEN
          FACT(0) = 1.D0