n  = 
sw  =  0.0;
swx  =  0.0;
swy  =  0.0;
swz  =  0.0;
varw  =  0.0;
varx  =  0.0;
vary  =  0.0;
varz  =  0.0;
idum  = 
ss  =  0.2*(exp(5.0)-exp(-5.0));
vol  =  3.0*7.0*ss;
for (j = 1; j<=n; j++)
{
    x = 1.0+3.0*ran2(idum);
    y = -3.0+7.0*ran2(idum);
    s = 0.00135+ss*ran2(idum);
    z = 0.2*log(5.0*s);
    if (z*z+pow(sqrt(x*x+y*y)-3.0,2) < 1.0)
	{
        sw = sw+1.;
        swx = swx+x;
        swy = swy+y;
        swz = swz+z;
        varw = varw+1.;
        varx = varx+x*x;
        vary = vary+y*y;
        varz = varz+z*z;
	}
}
w = vol*sw/n;
x = vol*swx/n;
y = vol*swy/n;
z = vol*swz/n;
dw = vol*sqrt((varw/n-power(sw/n,2))/n);
dx = vol*sqrt((varx/n-power(swx/n,2))/n);
dy = vol*sqrt((vary/n-power(swy/n,2))/n);
dz = vol*sqrt((varz/n-power(swz/n,2))/n);