function I = nestedintegration()

I = quad(@ghat, -1, 1) 


function z = ghat(y)

% y is a vector of evaluation points from quad.

z = zeros(size(y));

for i = 1:length(y),
    ybar = y(i);

% Contribution to the integral is zero if (1-ybar^2) <= 0.

    if ((1-ybar^2)>0)
        z(i) = quad(@(x)h(x,ybar), -sqrt(1-ybar^2), sqrt(1-ybar^2));
    end   
end   


function h = h(x,y)

hsq = 1-x.^2-y^2;

% Set contributions from outside the domain to zero.

hsq = (hsq > 0).*hsq; 
h = sqrt(hsq);