//==================================================== file = poisson.c =====
//=  Program to solve for Poisson probability                               =
//===========================================================================
//=  Notes: 1) Input from command line and output to stdout                 =
//=             * Input is three values, arrival rate lambda, service       =
//=               rate mu, and the number of events                         =
//=             * Ouput is Poisson probablity for number of events          =
//=         2) Uses an interative method to solve for Poisson probability   =
//=            for numerical stability                                      =
//=-------------------------------------------------------------------------=
//=  Example execution: (command line is poisson 600.0 100.0 10)            =
//=                                                                         =
//=    =============================================== poisson.c =====      =
//=    = lambda = 600.000000                                                =
//=    = mu     = 100.000000                                                =
//=    = m      = 10                                                        =
//=    ---------------------------------------------------------------      =
//=    = rho    = 6.000000                                                  =
//=    ---------------------------------------------------------------      =
//=    = Prob[number events <= m] = 0.957379                                =
//=    = Prob[number events == m] = 0.041303                                =
//=    = Prob[number events > m]  = 0.042621                                =
//=    ===============================================================      =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 poisson.c, cl poisson.c, gcc poisson.c                    =
//=-------------------------------------------------------------------------=
//=  Execute: poisson lambda mu m                                           =
//=-------------------------------------------------------------------------=
//=  Author: Ken Christensen                                                =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (08/26/17) - Genesis (from erlang.c)                      =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>          // Needed for printf()
#include <stdlib.h>         // Needed for exit(), atoi(), and atof()
#include <math.h>           // Needed for exp()

//===== Main program ========================================================
void main(int argc, char *argv[])
{
  double     lambda;        // Total arrival rate
  double     mu;            // Service rate
  int        m;             // Number of events to solve for
  double     rho;           // Utilization
  double     prob1;         // Prob[number events <= m]
  double     prob2;         // Prob[number events == m]
  double     prob3;         // Prob[number events > m]
  double     sum;           // Running sum value
  double     new_value;     // New value computed
  double     last_value;    // Last value computed
  int        i;             // Loop index

  // Check if sufficient command line parameters
  if (argc < 4)
  {
    printf("Usage is 'poisson lambda mu m' where lambda is the overall \n");
    printf("arrival rate, mu the service rate, and m is the number of  \n");
    printf("events to solve for.                                       \n");
    exit(1);
  }

  // Assign parameters from entered argv[]
  lambda = atof(argv[1]);
  mu = atof(argv[2]);
  m = atoi(argv[3]);

  // Compute probailities
  rho = lambda / mu;
  sum = last_value = 1.0;
  for (i=1; i<=m; i++)
  {
    new_value = last_value * (rho / i);
    sum = sum + new_value;
    last_value = new_value;
  }
  prob1 = exp(-rho) * sum;
  prob2 = exp(-rho) * last_value;
  prob3 = 1.0 - prob1;

  // Output input parameters and computed probabilities
  printf("=============================================== poisson.c =====\n");
  printf("= lambda = %f \n", lambda);
  printf("= mu     = %f \n", mu);
  printf("= m      = %d \n", m);
  printf("---------------------------------------------------------------\n");
  printf("= rho    = %f \n", rho);
  printf("---------------------------------------------------------------\n");
  printf("= Prob[number events <= m] = %f \n", prob1);
  printf("= Prob[number events == m] = %f \n", prob2);
  printf("= Prob[number events > m]  = %f \n", prob3);
  printf("===============================================================\n");
}