Numerically calculate the polar integral over angles I would like to calculate the integral of a data set defined for all angles (ie $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$). This data set is power in a given direction (power per unit solid angle) $P(\theta,\phi)$. The data set looks like (no, this question is not about parsing data files):
Theta Phi Power (W)
0     0   1.0
0     1   1.0
...
0     359 1.0
1     0   1.5
1     1   1.5
...
89    359 999.5
90    0   1000.0
90    1   999.5
...
180   359 1.0

What I have been doing is (knowing that $\Delta\theta=\Delta\phi=1\text{deg}=\frac{\pi}{180}\text{rad}$) calculating the integral
$$P_\text{tot}=\int_0^{2\pi}\!\int_0^\pi\! P(\theta,\phi) \sin \theta\; \text{d}\theta\; \text{d}\phi\approx\Delta\theta\Delta\phi\sum_{j=0}^{359}\sum_{k=0}^{180} P(k,j)\sin \left(k\frac{\pi}{180}\right)$$
But I am not certain that this is correct.
 A: The comments to the question have been helpful. The approximation is (almost) correct: I need to change the limit of my $\theta$ approximation sum to $k=179$:
$$  P_\text{tot}\approx\Delta\theta\Delta\phi\sum_{j=0}^{359}\sum_{k=0}^{179}P_{k,j}\sin\left(\frac{\pi k}{180}\right)$$
A: While $P(\theta,\,\phi)$ is from data in this case, it is also possible to use Monte Carlo methods for integrating; the usefulness of which is through it's constant order of error of $1/\sqrt{N}$, regardless of dimension. In general, the integral of $f(x)$ can be approximated as,
$$
I=\int f(x)\,\mathrm dx\approx\frac{1}{N}\sum_{i=1}^N\frac{f(x_i)}{p(x_i)}
$$
where $p(x)$ is a probability such that $\int_a^b p(x)\,\mathrm dx=1$ (e.g., in a uniform case, $p(x)=1/(b-a)$). Typically one tries to choose $p(x)\propto\vert f(x)\vert$, but that is not always possible.
In this case, your function would be something along the lines of 
function mc_integral(nmax):
    sum = 0.0
    for i =  1 to nmax do
        theta = arccos(1 - 2 * rand())
        phi = 2*pi*rand()
        sum = sum + P(theta, phi) * sin(theta)
    loop
    return sum / nmax * normalization

Through experimentation, it seems normalization=16.0, though at the moment I can't seem to figure out why. 
However, for your problem, since you have analytic data, you could amend the call to P(theta,phi) by converting those radians to (integer) angles and using a look-up table (2D array), p[theta_deg, phi_deg]. This may be more complicated than simply using the double sum as described in your post & answer, but it is still an alternative.
