src/math/RealSphericalHarmonic.cpp

#include <stdio.h>
#include <cmath>

#include "math/RealSphericalHarmonic.hpp"

using namespace oopse;

RealSphericalHarmonic::RealSphericalHarmonic() {
}

double RealSphericalHarmonic::getValueAt(double costheta, double phi) {
  
  double p, phase;
  
  // associated Legendre polynomial
  p = LegendreP(L,M,costheta);
  
  if (functionType == SH_SIN) {
    phase = sin((double)M * phi);
  } else {
    phase = cos((double)M * phi);
  }
  
  return coefficient*p*phase;
  
}

//---------------------------------------------------------------------------//
//
// double LegendreP (int l, int m, double x);
//
// Computes the value of the associated Legendre polynomial P_lm (x)
// of order l at a given point.
//
// Input:
//   l  = degree of the polynomial  >= 0
//   m  = parameter satisfying 0 <= m <= l,
//   x  = point in which the computation is performed, range -1 <= x <= 1.
// Returns:
//   value of the polynomial in x
//
//---------------------------------------------------------------------------//
double RealSphericalHarmonic::LegendreP (int l, int m, double x) {
  // check parameters
  if (m < 0 || m > l || fabs(x) > 1.0) {
    printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x);
    return NAN;
  }
  
  double pmm = 1.0;
  if (m > 0) {
    double h = sqrt((1.0-x)*(1.0+x)),
      f = 1.0;
    for (int i = 1; i <= m; i++) {
      pmm *= -f * h;
      f += 2.0;
    }
  }
  if (l == m)
    return pmm;
  else {
    double pmmp1 = x * (2 * m + 1) * pmm;
    if (l == (m+1))
      return pmmp1;
    else {
      double pll = 0.0;
      for (int ll = m+2; ll <= l; ll++) {
        pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
        pmm = pmmp1;
        pmmp1 = pll;
      }
      return pll;
    }
  }
}

Revision:	98
Committed:	Mon Oct 18 16:30:04 2004 UTC (20 years, 6 months ago) by gezelter
File size:	1768 byte(s)
Log Message:	Your basic run of the mill real-valued versions of Y_l^m(\theta, \phi)
#	User	Rev	Content
1	gezelter	98	#include <stdio.h>
2			#include <cmath>
3
4			#include "math/RealSphericalHarmonic.hpp"
5
6			using namespace oopse;
7
8			RealSphericalHarmonic::RealSphericalHarmonic() {
9			}
10
11			double RealSphericalHarmonic::getValueAt(double costheta, double phi) {
12
13			double p, phase;
14
15			// associated Legendre polynomial
16			p = LegendreP(L,M,costheta);
17
18			if (functionType == SH_SIN) {
19			phase = sin((double)M * phi);
20			} else {
21			phase = cos((double)M * phi);
22			}
23
24			return coefficientpphase;
25
26			}
27
28			//---------------------------------------------------------------------------//
29			//
30			// double LegendreP (int l, int m, double x);
31			//
32			// Computes the value of the associated Legendre polynomial P_lm (x)
33			// of order l at a given point.
34			//
35			// Input:
36			// l = degree of the polynomial >= 0
37			// m = parameter satisfying 0 <= m <= l,
38			// x = point in which the computation is performed, range -1 <= x <= 1.
39			// Returns:
40			// value of the polynomial in x
41			//
42			//---------------------------------------------------------------------------//
43			double RealSphericalHarmonic::LegendreP (int l, int m, double x) {
44			// check parameters
45			if (m < 0 \|\| m > l \|\| fabs(x) > 1.0) {
46			printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x);
47			return NAN;
48			}
49
50			double pmm = 1.0;
51			if (m > 0) {
52			double h = sqrt((1.0-x)*(1.0+x)),
53			f = 1.0;
54			for (int i = 1; i <= m; i++) {
55			pmm = -f h;
56			f += 2.0;
57			}
58			}
59			if (l == m)
60			return pmm;
61			else {
62			double pmmp1 = x * (2 * m + 1) * pmm;
63			if (l == (m+1))
64			return pmmp1;
65			else {
66			double pll = 0.0;
67			for (int ll = m+2; ll <= l; ll++) {
68			pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
69			pmm = pmmp1;
70			pmmp1 = pll;
71			}
72			return pll;
73			}
74			}
75			}
76