polar interpolation of data acquired in a cartesian plane

I'm working on a problem that I've modeled as:

[x,y] = meshgrid(-20:0.1:20,-20:2:20);
xx = vectorizeMatrix(x);
yy = vectorizeMatrix(y);
zz = zeros(size(xx));

r = 13;
A = r^2*pi;
zz(r^2>=xx.^2+yy.^2) = 1;
F = TriScatteredInterp(xx,yy,zz);
interpFun = @(xi,yi)F(xi,yi);
Aquad = quad2d(interpFun,-20,20,-20,20);
diff = percentDiff(A,Aquad,'type','first');

where 'vectorizeMatrix' and 'percentDiff' are functions of mine that, well, vectorize a matrix and calculate a percent difference.

I'm trying to very accurately calculate the area of the circle with the data using quad2d. The edges of the circle are giving quad2d problems, I believe, because of how it is interpolating the data (TriScatteredInterp, linear). It does relatively well at calculating the area, in this case it's <1% off, but I'd like it to get even better than that. Does anyone have any ideas about how to perform a polar interpolation on this data?