      PROGRAM PROBAB
************************************************************************
* Program to calculate the Poisson probability function using the      *
* function subprogram POISS.  Identifiers used are:                    *
*     LAMBDA : average # of occurrences of phenomenon per time period  *
*     N      : number of occurrences in a time period                  *
*     PROB   : Poisson probability                                     *
*     POISS  : function subprogram to calculate Poisson probability    *
*                                                                      *
* Input:   LAMBDA and N                                                *
* Output:  The Poisson probability PROB                                *
************************************************************************

      REAL LAMBDA, POISS, PROB
      INTEGER N

      PRINT *, 'THIS PROGRAM CALCULATES THE POISSON PROBABILITY FOR'
      PRINT *, 'LAMBDA = AVERAGE # OF OCCURRENCES PER TIME PERIOD'
      PRINT *, 'N = # OF OCCURRENCES FOR WHICH PROBABILITY TO BE FOUND'
      PRINT *, 'ENTER LAMBDA AND N (NEGATIVE VALUES TO STOP)'
      READ *, LAMBDA, N

* While LAMBDA >= 0 do the following

 10   IF (LAMBDA .GE. 0) THEN
         PROB = POISS(LAMBDA, N)
         PRINT 20, PROB
 20      FORMAT (1X, 'POISSON PROBABILITY = ', F6.4 /)
         PRINT *, 'ENTER LAMBDA AND N (NEGATIVE VALUES TO STOP)'
         READ *, LAMBDA, N
      GO TO 10
      END IF

      END

**POISS*****************************************************************
* Function to calculate the Poisson probability                        *
*                            N   -LAMBDA                               *
*                      LAMBDA * e                                      *
*          POISS(N) =  ------------------                              *
*                             N!                                       *
* Function FACTOR is called to calculate N!                            *
*                                                                      *
* Accepts:  Real number LAMBDA and integer N                           *
* Returns:  The poisson probability given by the formula above         * 
************************************************************************

      FUNCTION POISS(LAMBDA, N)

      INTEGER N, FACTOR
      REAL POISS, LAMBDA

      POISS = (LAMBDA ** N * EXP(-LAMBDA)) / FACTOR(N)
      END

**FACTOR****************************************************************
* Function to calculate the factorial N! of N which is 1 if N = 0,     *
* 1 * 2 * ... * N for N > 0.  It uses variable I as a  counter.        *
*                                                                      *
* Accepts:  Integer N                                                  *
* Returns:  The integer N!                                             *
************************************************************************

      FUNCTION FACTOR(N)

      INTEGER FACTOR, N, I

      FACTOR = 1
      DO 10 I = 2, N
         FACTOR = FACTOR * I
 10   CONTINUE
      END
