# Math Help - Numerical integration over a hemisphere

1. ## Numerical integration over a hemisphere

Dear all,

I am trying to solve the following problem:
A function is given on the surface of a unity hemisphere, and I need to find an efficient way to numerically integrate it over the surface.

I started with an equidistant approach (both angles of the spherical coordinates subdivided into constant intervals). This gives good results but needs a large amount of computations.
In the internet, I found an alternative (Gaussian) approach with equidistant points, but only for integration over spheres, not over hemispheres. This approach can be found here: Integration nodes for the sphere

Can somebody help me on doing the numerical integration over the hemisphere? Any help is appreciated!

Thank you,
Jens

2. Not sure what you want but if you choose to use Mathematica, here's the code I'd use to numerically integrate the function $f(x,y,z)$ over the upper hemisphere of the unit sphere.

Code:
NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1},
{y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]
I realize not everyone has Mathematica or even wants to use it. Just sayin' if you're interested in a brute force approach to the answer using Mathematica, which I'd think would be rather quick (at 2GHz) for a reasonable function, that's what I'd start with. Here's a check of $f(x,y,z)=x^2$ over the upper hemisphere of the unit sphere, the answer of which is $2\pi/3$, just to get some indication if the coding is correct:

Code:
f[x_, y_, z_] := x^2;
NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1},
{y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]
N[2*(Pi/3)]

2.0944

2.0944
That took lest than 1/10 of a second. Can you post your function?

3. shawsend,

I actually can't use Matematica in this case, because I need to include the integration (even for one special function only) into my own code. I'd like to find some algorithm and then implement it into my program for this one special case only.

Well, here is the function:
you are given 4 angles (describing 2 points on a unit hemisphere):
$\gamma, \varphi \in [-\pi, \pi]$
$\beta, \theta \in [0,\pi]$
$\gamma$ and $\theta$ are considered variables, while $\varphi$ and $\beta$ are given parameters.

The function I would like to integrate is the cosine of the angle between the position vectors of the 2 points:
$f_1(\gamma, \theta) = \cos\theta \sin\beta \cos(\gamma-\varphi) + \sin\theta \cos\beta$
Let $f_2(\gamma, \theta) = max[0,f_1(\gamma, \theta)]$, meaning that all angles $>\pi/2$ are discarded.

$\theta$ is the solar altitude angle. Since it is mathematically more common to start counting the angle from the z-axis, not from the x-y-plane, the following transformation is performed: $\theta_z = \pi/2-\theta$. This transforms the solar altitude angle $\theta$ into the zenith angle $\theta_z$. All following calculations are based on the zenith angle. In terms of $\theta_z$, the functions $f_1$ and $f_2$ read:
$f_1(\gamma,\theta_z) = \sin\theta_z \sin\beta \cos(\gamma-\varphi) + \cos\theta_z \cos\beta$
$f_2(\gamma, \theta_z) = max[0,f_1(\gamma, \theta_z)]$

Then, the integration I am trying to calculate numerically is:
$\phi = g(\beta,\varphi) = \frac{1}{\pi}\int_{\gamma=-\pi}^{\pi}{\int_{\theta_z=0}^{\pi/2}{f_2(\gamma, \theta_z)}\sin\theta_z\ d\theta_z d\gamma}$
I hope that you can now better see what I mean and what I want to do.
Jens

4. Ok, I see what you mean. I've done numerical integration in C++. Just easier with Mathematica. Also, can link Mathematica to an external program but that's a programming project and something you're probably not interested in. I may still work with your integral to see what I get with Mathematica.

5. thanks for your effort!
i am looking forward to getting your results.

Jens

6. Bamus, your notation is unusual to me (not the standard way to set up spherical coordinates and dot product of two vectors, at least to me). In your notation, the point $p(\gamma,\beta)$, $\gamma$ is measured in the counter-clockwise direction from the negative x-axis, and $\beta$ is measured up from the x-y plane? Same with the other two symbols?

Also, your integral looks a bit odd as your integrating $f_2(\gamma,\theta)$ but using $d\theta d\phi$. I would gusss perhaps it would be something like:

$g(\beta,\varphi)=\frac{1}{\pi}\int_{-\pi}^{\pi}\int_{0}^{\pi} f_2(\gamma,\theta)sin(\theta)d\theta d\gamma$

Maybe no though. Anyway, generating a plot for $g(\beta,\varphi)$ is problematic as the integration is taking a long time to complete. Probably due to the complicated integrand. So I switched to Monte Carlo method. This was better. Perhaps you might want to try this method in your approach to the problem

7. shawsend,

thanks again.
Sorry for the mistakes in the formulae and thank you for pointing them out. I corrected the errors, see my edited post above.

Actually, I do not need to have the value of $\phi$, the integral, calculated as a general function $\phi(\beta, \varphi)$. $\beta$ and $\varphi$ are just parameters, not variables in the proper sense. I am given values for $\beta$ and $\varphi$ at the start of the calculation, and I am supposed to do the calculation of $\phi$ only once and only for this pair of $\beta$ and $\varphi$.

So, for the purpose of calculating the integral, we can treat $\beta$ and $\varphi$ as if they were real constants. Will we still need Monte Carlo? I have never worked with that method...

Jens

8. Hello Bamus,

How about posting your values of $\beta,\varphi$. I'd like to generate a more accurate integral and see if it agrees with your value. Also Monte Carlo is a random method. Just like if you wanted to calculate the area under a parabola in the unit square, generate 100,000 random points in that square and count the number inside the parabola. The ratio inside/outside gives a pretty accurate value to the integral. The method is commonly used to calculate many-iterated integrals of many variables.

That's just one method out of many available to Mathematica, but I would think that's unnecessary since we're only calculating one point and not 10,000 or so to make a decent plot of $g(\beta,\varphi)$.

9. -> function $f_1(\gamma, \theta)$ corrected and simplified a little bit, see post above.

Originally Posted by shawsend
How about posting your values of $\beta,\varphi$. I'd like to generate a more accurate integral and see if it agrees with your value.
For instance, typical values for $\beta$ and $\varphi$ would be:
$\beta = 35$° $= 0.61086524$ rad
$\varphi = 15$° $= 0.26179939$ rad

You could check with
$\beta = 0$° and $\varphi = 0$° which must give $g(0,0) = 1$
$\beta = 90$° and $\varphi = 0$° which must give $g(\pi/2,0) = 1/2$

Originally Posted by shawsend
That's just one method out of many available to Mathematica, but I would think that's unnecessary since we're only calculating one point and not 10,000 or so to make a decent plot of $g(\beta,\varphi)$.
That's why in a recent post I pointed out that I don't need a general function $g(\beta,\varphi)$.

Thanks for your help!
Jens

10. Then I must be doing something wrong because I don't get those answers for $g(0,0]$ and $g(\pi/2,0)$. The code below is a little messy because of the greek letters but you can clearly see the last two command are the results for these two values.

f[\[Gamma]_, \[Beta]_, \[CurlyPhi]_, \[Theta]_] :=
Max[0, Cos[\[Theta]]*Cos[\[Gamma] - \[CurlyPhi]]*Sin[\[Beta]] +
Cos[\[Beta]]*Sin[\[Theta]]]
g[\[Beta]_, \[CurlyPhi]_] := (1/Pi)*NIntegrate[f[\[Gamma], \[Beta], \[CurlyPhi], \[Theta]]*
Sin[\[Theta]], {\[Theta], 0, Pi}, {\[Gamma], -Pi, Pi},
Method -> MonteCarlo];
g[0, 0]
g[Pi/2, 0]

Out[3]=
3.1486450543682474

Out[4]=
0.6402232798879435

However I did generate a 15x15 array of points for $g(\beta,\varphi)$ and plotted them. Is the plot below consistent with your expectations?

11. I apologize again... I have been confusing about the use of the solar altitude angle $\theta$ vs. the zenith angle $\theta_z$. This changes some $\sin$ into $\cos$ functions (and vice versa) in the functions. See the details and the corrected function in my post above, containing the definition of $f_1(\gamma, \theta_z)$.

The two tests yield:
Special case 1:
$\beta = 0$°, $\varphi = 0$°:
In this case, $f_1(\gamma,\theta_z)=\cos \theta_z$, and as $f_1\ge0$ is always true for $\theta_z \in [0,\pi/2]$, the value of the integral is: $g(0,0) = 1$ (checked this by hand).

Special case 2:
$\beta = 90$°, $\varphi = 0$°:
In this case, $f_1(\gamma,\theta_z)=\sin \theta_z \cos\gamma$. Here, the condition $f_1(\gamma,\theta_z)\ge0$ is equivalent to $\gamma \in [-\pi/2,\pi/2]$. Thus, the value of the integral $g(\pi/2,0) = 1/2$ can be easily calculated by hand.
Sorry again for the confusion.
Now the problem should be well-defined...

Jens