How can i evaluate this surface integral? (It has some singularities)

Im trying to integrate the following function over X (0->Pi) and then y (0-> 2*Pi).

When i create a surface (x,y,z) after a very fine meshgrid (1000 x 1000), i see what appears to be 4 very sharp singularities.

Is there any way to evaluate this integral by ignoring or smoothing over those singularities?

Any help would be greatly appreciated.

FYI i am using matlab 2011b.

@(x,y)-sin(x).*(sin(y).^2.*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(cos(x).^2.*cos(y).^2.*1.05e2+cos(y).^2.*sin(x).^2.*4.58e2-cos(x).^4.*sin(y).^2.*9.72e2-sin(x).^4.*sin(y).^2.*9.72e2+cos(x).^2.*sin(x).^2.*sin(y).^2.*1.944e3).*(3.0./5.0e3)+cos(x).^2.*cos(y).^2.*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(cos(x).^2.*cos(y).^2.*-1.087e3+cos(y).^2.*sin(x).^2.*3.15e2+cos(x).^4.*sin(y).^2.*3.15e2+sin(x).^4.*sin(y).^2.*3.15e2-cos(x).^2.*sin(x).^2.*sin(y).^2.*6.3e2).*(1.0./5.0e3)+cos(y).^2.*sin(x).^2.*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(cos(x).^2.*cos(y).^2.*1.05e2-cos(y).^2.*sin(x).^2.*9.72e2+cos(x).^4.*sin(y).^2.*4.58e2+sin(x).^4.*sin(y).^2.*4.58e2-cos(x).^2.*sin(x).^2.*sin(y).^2.*9.16e2).*(3.0./5.0e3)-cos(x).^2.*cos(y).^4.*sin(x).^2.*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(4.29e2./5.0e2)-cos(x).*cos(y).^2.*sin(y).*(cos(x).^3.*sin(y)-cos(x).*sin(x).^2.*sin(y)).*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(4.29e2./5.0e2)+cos(y).^2.*sin(x).*sin(y).*(sin(x).^3.*sin(y)-cos(x).^2.*sin(x).*sin(y)).*1.0./(cos(x).^2.*cos(y).^2+cos(x).^2.*sin(y).^2+cos(y).^2.*sin(x).^2-sin(x).^2.*sin(y).^2).^2.*(4.29e2./5.0e2))

Answer 1

[Mike Hosea]

When I split the integral up into regions with difficult parts (a couple of circles) on the boundaries, I can get QUAD2D and INTEGRAL2 to integrate over the regions only if I cap the values of the integrand function (z(z>M) = M and z(z<-M)=-M). Because the value of the integral increases steadily as I increase M (until the integration fails because of numerical or minimum step size issues), I don't think the singularities are integrable, i.e. whatever you do to mitigate the singularities will probably just change the problem to an easier problem with a different answer.