Solid Angle Integrals

Question

I'm currently working on solving a solid angle calculation for my physics project. I'm trying to solve it computationally using python - but my coding skills are basic at best. I need some hints/tips on how to go about solving this equation which I obtained from a research paper - equation 9b on (http://ac.els-cdn.com/S0969804306002090/1-s2.0-S0969804306002090-main.pdf?_tid=7b61d61e-d162-11e4-9032-00000aab0f6c&acdnat=1427118343_91eaa37dbd9242927e86c41952074750)

The integral is:

$$\int_{12.76}^{43.71}\text{arccos}\left(\frac{p^2 - r^2 + h^2\tan^2(x)}{2ph\tan(x)}\right)\sin(x)\; dx$$

with $h = 68,$ $r = 24.8$ and $p = 40.2$

Any help appreciated, thanks.

Answer 1

For this answer I assume that everything is in degrees, though of course that's just a matter of scaling.

When you want to do numerical computing in python there's two handy packages: numpy and scipy. Using these, your computation takes only a few lines. First define a function for your integrand:

def integrand(x):
  import numpy as np
  p = 40.2
  r = 24.8
  h = 68
  theta = np.pi * x / 180 
  numerator = p**2 - r**2 + (h**2) * np.tan(theta)**2
  denominator = 2*h*p*np.tan(theta)
  return (np.arccos(numerator/denominator) * 180/np.pi) * np.sin(theta)

This could have been written in fewer lines, I just found that this way was easier to debug. Now to integrate, you simply run

from scipy import integrate
result = integrate.quad(integrand, 12.77, 43.7)
print(result)

I had to truncate the interval a tiny bit because at the endpoints the input to $\arccos$ evaluates to slightly over $1$ and this ruins the whole computation.

The answer I got was about $425.4391$, to about that much accuracy (though this is effected by the truncation I did, so don't trust that number too much!).