# 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.

• Looks OK, but what about trying out some non-trivial $P$'s, where you know the answer? Feb 7, 2017 at 9:22
• In your data file you seem to have $\theta$ running from 0 to 359, not to 360 as in your approximation. Also, I believe that if you divide the interval in $N$ pieces you should have a sum of $N$ values, i.e. if you divide the integration over $\phi$ in 180 pieces, you should have the sum over 180 (not 181) values. Feb 7, 2017 at 9:42
• Being nitpicky, but you really want $\sin(k\frac{\pi}{180})$, unless you mean a version of $\sin$ which takes in degrees instead of radians. Make sure that if $P(\theta,\phi)=1$ that you get very close to $4 \pi$. Feb 7, 2017 at 9:45
• @all Thanks for your comments. I have run this for $P(\theta,\phi)=$ and $P(\theta,\phi)=\sin\theta|\cos\phi|$ and have respectively got $4\pi$ and $2\pi$, indicating that this approximation is working. And yes, I should really be summing up to $k=179$ on the inner sum as I've split the $\theta$ integral up into 180 pieces. Feb 7, 2017 at 23:35

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)$$

• Very nice. Might also be useful to look into Monte Carlo techniques for integrating over solid angles. Feb 8, 2017 at 2:46
• Right, that would be useful given the $P(\theta,\phi)$ as a function, but I have $P_{k,j}$ as a data set, and I need to integrate over the entire solid angle range and get the total (power) Feb 8, 2017 at 12:26
• So you're saying that won't work? I might disagree with that. Feb 8, 2017 at 13:06
• Nah, not saying it won't work, just not worth the effort given I have a simple data set defined for 181 by 361 values. Though I'm not sure how you would implement Monte-Carlo methods on a data set, as I've only learnt how to use those methods given the mathematical function expression Feb 8, 2017 at 13:08
• Well that's a more reasonable criticism ;) Feb 8, 2017 at 13:38

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.