double poidev(int xm, int &idum)
{
	double pi,g,em,t,sq,alxm,y,dum,temp;
	int oldm=0;
    pi = 3.141592654;
    if (xm < 12) 
	{
        if (xm != oldm)
		{
            oldm = xm;
            g = exp(-xm);
        }
        em = -1.0;
        t = 1.0;
        do
		{
			em = em + 1.0;
			t = t * ran1(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 * ran1(idum));
				em = sq * y + xm;
			}
			while (em < 0.0);
			em = int(em);
			dum = em * alxm - gammln(em + 1.0) - g;
			t = 0.9 * (1.0 + y *y) * exp(dum);
		}
        while (ran1(idum) > t);
    }
    temp = em;
	return temp;
}