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$
  • $\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 Mar 30 '15 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 Mar 30 '15 at 14:27
  • $\begingroup$ I checked different spherical harmonics and they seem fine - only $Y^{0}_{2}$ is bad! $\endgroup$ – WoofDoggy Mar 30 '15 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 Mar 31 '15 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 Mar 31 '15 at 15:21
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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.