/********************************************************************/
/* 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);
}