function p = walk_on_spheres(a,b,m,h)

n = 1e6; N = m*n;
p = 0;
for k=1:m
    x = zeros(1,n); y = x;
    while ~isempty(x)
        dx = min(a-x,a+x); dy = min(b-y,b+y);
        r = min(dx,dy);
        phi = 2*pi*rand(size(x));
        x = x+r.*cos(phi); y = y+r.*sin(phi);
        tx = abs(x)>a-h; ty = abs(y)>b-h;
        ts = ~(tx|ty);
        p = p + sum(tx);
        x = x(ts); y=y(ts);
    end
end
p = p/N;
return;