double gamdev(int ia, int &idum)
{
	double x,v1,v2,y,s,e,t;
	int j,am;
    if (ia < 1)
	{
		cout<< "abnormal exit";
		exit(1);
	}
    if (ia < 6)
	{
        x = 1.0;
        for (j = 1; j<=ia; j++)
            x = x * ran1(idum);
        x = -log(x);
	}
    else
	{
		do
		{
			do
			{
			do
			{
				v1 = 2.0 * ran1(idum) - 1.0;
				v2 = 2.0 * ran1(idum) - 1.0;
			}
			while (v1*v1 + v2*v2 > 1.0);
			y = v2 / v1;
			am = ia - 1;
			s = sqrt(2.0* am + 1.0);
			x = s * y + am;
			}
			while (x <= 0.0);
			e = (1.0 + y *y) * exp(am * log(x / am) - s * y);
		}
		while (ran1(idum) > e);
	}
    t = x;
	return t;
}