double gasdev(int& idum)
{
	static int iset;
	static double gset;
	double t,v1,v2,fac,r;
    if (iset == 0.0)
	{
		do
		{
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			r = v1 *v1 + v2 *v2;
		}
		while ((r >= 1.0) || (r == 0.0));
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		t = v2 * fac;
		return t;
		iset = 1;
	}
    else
	{
		t = gset;
		return t;
		iset = 0;
    }
}