Numerical integration over 2D disk I have a real-valued function $f$ on the unit disk $D$ that is fairly well behaved (real-analytic everywhere) and would like to find the integral $\int_D f(x,y)dxdy$ numerically. After much searching, the best methods I've been able to find have been the simple quadrature rules in Abramowitz and Stegun that sample $f$ at up to 21 points. What work has been done since? In particular, I'm interested in rules that allow sampling at more than 21 points. One reference indicates that finding optimal quadrature rules is a hard problem, but it seems to me that something better must have been published in the last 50 years.
A couple of references suggest integrating over various domains by triangulating them and using numerical integrals over the triangles. Is this the preferred method for a disk?
(I'm trying to improve code that has been implemented using quasi-monte carlo. It seems to me that we could do much better using the knowledge that $f$ is real-analytic and probably well approximated by polynomials.)
Update: I can't easily say exactly what the functions $f$ are as they're the messy result of a chain of computations. I can say that qualitatively it's like a gaussian with a central hump, fast decay, though not exact rotational symmetry. I do have a pretty good handle on how big the hump is and where it's centred. All variations on this might happen and some are easy to dispense with: eg. the hump might be situated well outside the disk so I know the integral is nearly zero. Or the hump may be very wide in which case the integrand is almost constant. Sometimes the hump is contained well within the disk in which case I can switch to more efficient quadrature over the (approximate) support of $f$ rather than the disk. But having said all that, I'd still like to see some general gaussian quadratures rules for the disk that would apply to integrating any function over the disk that is well approximated by a polynomial.
Update2: After much web searching I found some Fortran code to do what I want (and more) and a reference to a book by Arthur Stroud, Approximate Calculation of Multiple Integrals. It seems as that this work from 1971 is the state of the art.
 A: See the Encyclopedia of Cubature Formulas.  The site is password protected, but the maintainer will give a password to anyone who asks.  
A: Here's another one. a little strange, since the quadrature rule involve the value of line integrals, but I guess you can write those using a separate quadrature rule...
Bojanov, Borislav; Petrova, Guergana, Numerical integration over a disc. A new Gaussian quadrature formula, Numer. Math. 80, No. 1, 39–59 (1998). ZBL0911.65015.

We construct a quadrature formula for integration on the unit disc which is based on line integrals over $n$ distinct chords in the disc and integrates exactly all polynomials in two variables of total degree $2n - 1$.

Cheers.
A: You might like this one, I do not know if the entire link will fit.
Wait, it is available for download from her website!
http://www.math.tamu.edu/~gpetrova/
Journal of Approximation Theory
 Volume 104, Issue 1,  (May 2000)
Uniqueness of the Gaussian Quadrature for a Ball
Pages 21-44
Borislav Bojanov and Guergana Petrova
Department of Mathematics, University of Sofia, Boulevard James Boucher 5, 1164, Sofia, Bulgariaf1
Department of Mathematics, University of South Carolina, Columbia, South Carolina, 29208, U.S.A., f2
Received 8 June 1999; 
accepted 22 October 1999. ;
Available online 26 March 2002.
Abstract
We construct a formula for numerical integration of functions over the unit ball in Image d that uses n Radon projections of these functions and is exact for all algebraic polynomials in Image d of degree 2n−1. This is the highest algebraic degree of precision that could be achieved by an n term integration rule of this kind. We prove the uniqueness of this quadrature. In particular, we present a quadrature formula for a disk that is based on line integrals over n chords and integrates exactly all bivariate polynomials of degree 2n−1.
Author Keywords: Gauss quadrature formula; orthogonal polynomials; highest degree of precision
Different article by same people:
http://www.math.tamu.edu/~gpetrova/CAM7238.pdf
A: Here is link to my site, where I publish some of the cubatures of such kind:
Cubature formulas for the unit disk
Page also gives C source code for the numerical integration by the product of two 1D Gauss-type quadratures. Images of points distribution in that case are available too.
Let me know if you find it useful.
A: An insufficiently well-known (so, perhaps slightly beyond the state of the art) integration algorithm can be found in the paper of O. Jenkinson and M. Pollicott entitled
"A dynamical approach to accelerating numerical integration with equidistributed points".
They claim a qualitative improvement over the state of the art as of three years ago.
A: An interesting approach can be found in the paper Extensions of Gauss Quadrature Via Linear Programming by Ernest Ryu and Stephen Boyd.
Their method is to show that ordinary Gauss quadrature on the real line can be interpreted as a solution to a linear programming (LP) problem and then use the same LP method on domains other than the line. When used in a more general setting it's no longer solving precisely the same problem that Gauss quadrature is solving but it's still a good heuristic for generating efficient quadrature methods.
One nice thing about their approach is that it generalizes to irregular domains too.
A: Forgive me for stating the obvious, but your best bet will be to use matlab/octave/scilab/.... 
I'd assume that any modern general purpose code will be adaptive (subdivide the area into squares recursively, estimate the integral on each square using a simple quadrature rule and
spend most effort on the areas that are converging most slowly). 
I only know enough to be dangerous. I would not attempt to write my own code for this purpose,
since the code produced by a real numerical analyst will be vastly superior.
