0
$\begingroup$

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:

1) Go to link and download source code .c and.h files

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

$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – flawr
    Commented Mar 30, 2015 at 13:20
  • $\begingroup$ 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) $\endgroup$
    – WoofDoggy
    Commented Mar 30, 2015 at 14:27
  • $\begingroup$ I checked different spherical harmonics and they seem fine - only $Y^{0}_{2}$ is bad! $\endgroup$
    – WoofDoggy
    Commented Mar 30, 2015 at 21:57
  • $\begingroup$ 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.) $\endgroup$
    – flawr
    Commented Mar 31, 2015 at 15:07
  • $\begingroup$ 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 :( $\endgroup$
    – WoofDoggy
    Commented Mar 31, 2015 at 15:21

1 Answer 1

0
$\begingroup$

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .