# Lebedev grid integral with $\cos(2\theta)$.

This is just a part of the problem I have noticed while dealing with numerical integration in spherical coordinates. I want to evaluate integral

$$\int\limits_{0}^{2\pi}d\varphi\int\limits_{0}^{\pi}d\theta\sin(\theta)\ (1 + 3\cos(2\theta))$$

using Lebedev grid. This integral is obviously $0$, but not in numerical scheme. I've implemeted this integral in C using lebedev library link with finest grid (5810), but it is not working. Not only this function, but just $\cos(2\theta)$. Basically for $\cos(2\theta)$ it gives me $-0.000019$ which is far from exact value $-0.66667*2\pi$.

In Lebediev method You replace integral of the form $$\int\limits_{0}^{2\pi}d\varphi\int\limits_{0}^{\pi}d\theta\sin(\theta)\ f(\theta,\phi)$$ with single sum $$4\pi\sum\limits_{i=0}^{N}w_{i}f(\theta_{i}, \phi_{i})$$ using weights $w_{i}$.

Observation:

I did a test on spherical harmonics $Y^{0}_{0}, Y^{0}_{1}, Y^{0}_{2}$ and first two are correct, but the third one is obviously wrong! As far as I know Lebedev rule should be exact for spherical harmonics.

Minimal working example:

2) Here is my program to integrate spherical harmonics $Y^{0}_{0}$ (Y00), $Y^{0}_{1}$ (Y10), $Y^{0}_{2}$ (Y20)

# include <stdlib.h>
# include <math.h>
# include <stdio.h>

# include "sphere_lebedev_rule.h"

# define SIZE 14  // size of the Lebedev grid
# define PI 3.14159265358979323846

int main () {

double *X, *Y, *Z, *Th, *Ph, *W;
double Y00 = 0.0, Y10 = 0.0, Y20 = 0.0;
int i;

// Allocate memeory for spherical grid
X = (double*) malloc (SIZE*sizeof(double));
Y = (double*) malloc (SIZE*sizeof(double));
Z = (double*) malloc (SIZE*sizeof(double));
Th = (double*) malloc (SIZE*sizeof(double));
Ph = (double*) malloc (SIZE*sizeof(double));
W = (double*) malloc (SIZE*sizeof(double));

// Generate grid
ld0014 (X, Y, Z, W);
// Prepare grid in spherical coordinates
for (i = 0; i<SIZE; i++) {
xyz_to_tp (X[i], Y[i], Z[i], Th+i, Ph+i);
}
// free cartesian coordinates
free (X);
free (Y);
free (Z);

// evaluate integrals of spherical harmonics
for (i = 0; i<SIZE; i++) {
Y00 += W[i];
Y10 += W[i]*cos(Th[i]);
Y20 += W[i]*(3.0*cos(Th[i])*cos(Th[i]) - 1.0);
}

Y00 *= 0.5*1.0/sqrt(PI)*4.0*PI;
Y10 *= 0.5*sqrt(3.0/PI)*4.0*PI;
Y20 *= 0.25*sqrt(5.0/PI)*4.0*PI;

printf ("Y00 = %.20f, Y10 = %.20f, Y20 = %.20f\n", Y00, Y10, Y20);

// free memory
free (Th);
free (Ph);
free (W);

return 0;
}


3) compile: gcc -Wall [nazwa].c sphere_lebedev_rule.c -lm

• Can you provide a minmal working example for us to try it ourselves? Have you also tried integrating simpler functions like $f(\theta, \phi) = c \in \mathbb R$ or so? – flawr Mar 30 '15 at 13:20
• I will put an example in C later, but simple constant function is ok, because You are summing over the weights ($Y_{0}^{0}$ example) – WoofDoggy Mar 30 '15 at 14:27
• I checked different spherical harmonics and they seem fine - only $Y^{0}_{2}$ is bad! – WoofDoggy Mar 30 '15 at 21:57
• I could not get it compiled as I have no experience with C. What happens for different grid sizes? Can you give a definition of $Y_k^n$? (The one I've seen I've never had anything to do with sphercial harmonics either.) – flawr Mar 31 '15 at 15:07
• Spherical harmonics are listed here. The larger the grid size the worse is Y20. It is really strange for me and I have no idea what I am doing wrong. I am considering taking different rule if this is not working :( – WoofDoggy Mar 31 '15 at 15:21

The problem lies in the library function:

void xyz_to_tp ( double x, double y, double z, double *t, double *p )

/******************************************************************************/
/*
Purpose:

XYZ_TO_TP converts (X,Y,Z) to (Theta,Phi) coordinates on the unit     sphere.

Modified:

09 September 2010

Author:

Dmitri Laikov

Parameters:

Input, double X, Y, Z, the Cartesian coordinates of a point
on the unit sphere.

Output, double T, P, the Theta and Phi coordinates of
the point.
*/
{
double ang_x;
double fact;
double pi = 3.14159265358979323846;

*p = acos ( z );

fact = sqrt ( x * x + y * y );

if ( 0 < fact )
{
ang_x = acos ( x / fact );
}
else
{
ang_x = acos ( x );
}

if ( y < 0 )
{
ang_x = - ang_x;
}
*t = ang_x;
/*
Convert to degrees.
*/
*t = *t * 180.0 / pi;
*p = *p * 180.0 / pi;

return;
}


In standard convention

$$x = \sin\theta \cos\varphi \\ y = \sin\theta \sin\varphi \\ z = \cos\theta$$

so,

$$\theta = \text{acos}(z) \\ \varphi = \text{acos}\left(\frac{x}{\sqrt{x^2 + y^2}} \right)$$ In library function $p$ is actually THETA and $t$ is PHI.