#include <iostream.h>
#include <math.h>
#include "ranpois.h"
#include "random.h"
double pois (int, double);

int ranpois (double nu, int& seed){

  /*  Author:  Glen Cowan
      Date:    30 Oct 2001
      Function to generate Poisson random numbers according to the
      algorithm in Sheldon M. Ross, A Course in Simulation, Macmillan, 
      1990, page 49.
  */

  int n_max = int ( nu * 100. );
  double r = random (seed);
  bool find_n = true;
  int n = 0;
  double v = pois(n, nu);
  double u = random (seed);
  while ( find_n ){
    if ( u < v ) {
      find_n = false;                                   // success
    } else if ( n > n_max ) {
      cout << "ranpois>> unable to generate n" << endl;
    } else {
           n = n + 1;
           v = v + pois (n, nu);
    }
  }

  return n;

}

double pois(int n, double nu){

//  returns Poisson probability for n given parameter nu.

  double p = exp ( - nu );
  for(int m = 1; m<=n; m++) p = p * nu / double(m);
  return p;
}