PROGRAM Poisson_Probability
!-----------------------------------------------------------------------
! Program to calculate the Poisson probability function using the
! function subprogram Poisson.  Identifiers used are:
!   AveOccurs   : average # of occurrences of phenomenon per
!                 time period
!   NumOccurs   : number of occurrences in a time period
!   Probability : Poisson probability
!   NumProbs    : number of probabilities to calculate
!   I           : DO-loop control variable
!   Poisson     : internal function subprogram to calculate Poisson
!                 probability
!   Factorial   : internal function subprogram to calculate factorials
!
! Input:   NumProbs and values for AveOccurs and NumOccurs
! Output:  Poisson probabilities
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: AveOccurs, Probability
  INTEGER :: NumProbs, I, NumOccurs

  PRINT *, "This program calculates Poisson probabilities."
  WRITE (*, '(1X, A)', ADVANCE = "NO") &
         "How many probabilities do you wish to calculate? "
  READ *, NumProbs

  DO I = 1, NumProbs
     WRITE (*, '(1X, A)', ADVANCE = "NO") &
           "Enter average # of occurrences per time period: "
     READ *, AveOccurs
     WRITE (*, '(1X, A)', ADVANCE = "NO") &
           "Enter # of occurrences for which to find probability: "
     READ *, NumOccurs
     Probability = Poisson(AveOccurs, NumOccurs)
     PRINT '(1X, "Poisson probability = ", F6.4 /)', Probability
  END DO

CONTAINS

  !-Poisson ------------------------------------------------------------
  ! Function to calculate the Poisson probability
  !                            N   -Lambda
  !                      Lambda * e
  !          Poisson(N) =  ------------------
  !                             N!
  ! Function Factorial is called to calculate N!
  !
  ! Accepts:  Real number Lambda and integer N
  ! Returns:  The poisson probability given by the formula above
  !---------------------------------------------------------------------
  
  FUNCTION Poisson(Lambda, N)
  
    REAL :: Poisson
    REAL, INTENT(IN) :: Lambda
    INTEGER, INTENT(IN) :: N
  
    Poisson = (Lambda ** N * EXP(-Lambda)) / REAL(Factorial(N))

  END FUNCTION Poisson
  
  
  !- Factorial ---------------------------------------------------------
  ! Function to calculate the factorial N! of N which is 1 if N = 0, 
  ! 1 * 2 * . . . * N if N > 0.
  !
  ! Accepts:  Integer N
  ! Returns:  The integer N!
  !
  ! Note:  I is a local integer variable used as a counter.
  !---------------------------------------------------------------------
  
  FUNCTION Factorial(N)
  
    INTEGER, INTENT(IN) ::  N
    INTEGER :: Factorial, I
  
    Factorial = 1
    DO I = 2, N
       Factorial = Factorial * I
    END DO

  END FUNCTION Factorial
  
END PROGRAM Poisson_Probability