Numerical integration over a hemisphere

**cool.bambus** · September 18th 2008, 04:36 AM

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

**shawsend** · September 18th 2008, 07:56 AM

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?

**cool.bambus** · September 18th 2008, 05:05 PM

shawsend,

Thank you for your answer!
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.
Thanks in advance for your help!
Jens

**shawsend** · September 18th 2008, 11:56 PM

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.

**cool.bambus** · September 19th 2008, 11:41 AM

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

Jens

**shawsend** · September 20th 2008, 05:46 AM

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

**cool.bambus** · September 20th 2008, 04:50 PM

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

**shawsend** · September 21st 2008, 04:28 PM

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

**cool.bambus** · September 21st 2008, 04:59 PM

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

**shawsend** · September 22nd 2008, 03:41 AM

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?

**cool.bambus** · September 22nd 2008, 07:53 AM

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

Thread: Numerical integration over a hemisphere

LinkBack

Thread Tools

Display

Numerical integration over a hemisphere

Similar Math Help Forum Discussions

numerical integration

Example of Numerical Integration?

Numerical integration

multiple integration: hemisphere/cylinder

Numerical Integration

Search Tags