/********************************************************************/
/*    rand.c                                                        */
/*                                                                  */
/*    Routine to return a random number from a Poisson distribution */
/*    Includes the following routines from Numerical Recipes in C,  */
/*    p 222-223, by Press WH, et al. Cambridge University Press (C) */
/*    1988:                                                         */
/*                                                                  */
/*    expdev(), poidev() and gammln()                               */
/*                                                                  */
/*    Please ensure you have the necessary licenses to use these.   */
/*                                                                  */
/********************************************************************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define PI      3.141592654

double gammln(double x);
double poidev(double xm, long *idum);
double expdev (double xm, long *idum);

double expdev(double xm, long *idum)
{
   double dum;

   do
      dum = ran2(idum);
   while (dum == 0.0);
   return (-xm*log(dum));
}

double poidev(double xm, long *idum)
{
 static double sq, alxm, g, oldm=(-1.0);
 double em, t, y;

 if (xm < 12.0)
 {
  if (xm != oldm)
  {
   oldm = xm;
   g = exp(-xm);
  }

  em = -1;
  t = 1.0;
  do {
   em += 1.0;
   t *= ran2(idum);
  } while (t > g);

 }
 else
 {
  if (xm != oldm)
  {
   oldm = xm;
   sq = sqrt(2.0*xm);
   alxm = log(xm);
   g = xm * alxm - gammln(xm+1.0);
  }
  do {
   do {
    y = tan(PI*ran2(idum));
    em = sq * y + xm;
   } while (em < 0.0);
   em = floor(em);
   t = 0.9 * (1.0 + y*y) * exp(em * alxm - gammln(em+1.0) - g);
  } while (ran2(idum) > t);

 }

 return em;

}


double gammln(double xx)
{
 double x, tmp, ser;
 static double cof[6]={76.18009173,-86.50532033,24.01409822,
	-1.231739516,0.120858003e-2,-0.536382e-5};
 int j;

 x = xx - 1.0;
 tmp = x + 5.5;
 tmp -= (x + 0.5) * log(tmp);
 ser = 1.0;
 for (j=0; j<=5; j++)
 {
  x += 1.0;
  ser += cof[j]/x;
 }

 return -tmp + log(2.50662827465 * ser);

}