I am struggling with a problem that requires me to numerical integrate a function within a hemisphere. While the function is smooth, it is pretty nasty, extremely oscillatory and I have to calculate it at any given point approximately from solving another problem in Stokes flow.
At the moment, I have just naively been integrating it in a hemisphere by using spherical coordinates and doing trapezoid in the $\theta$ dimension and Gaussian in the $p$ and $\phi$ dimensions. The numerical answer I get out doesn't settle down until roughly 500,000-1 million quadrature points are used (it's pretty similar if I'm integrating over a cube and using Gaussian quadrature in all the dimensions). So a few questions
1) If I'm doing spherical, then the points are tending to cluster around the origin with my naive approach. What is a good way to weight the node in the $p$ dimension so they are more evenly distributed? I did some Googling and saw one suggestion of using a weighted Gaussian with a weight $x$, but I don't know a priori why that would be the best choice, and I can't find any clear lists or programs for Matlab for weighted Gaussian nodes and if I'm writing it from scratch I might as well use the most 'ideal' weight function I can.
2) Is it a terrible idea to just throw more Gaussian points at it? Should I instead be cutting up the region into smaller and smaller sections and using the same number of quadrature points in each region? If so, I'm not sure how to subdivide my hemispherical region in a reasonable way.
3) Any other suggestions?