trunk/SHAPES/SHFunc.cpp

#include <stdio.h>
#include <cmath>
#include "SHFunc.hpp"

SHFunc::SHFunc() {
}

double SHFunc::getValueAt(double costheta, double phi) {

  double p, phase;

  // associated Legendre polynomial
  p = LegendreP(L,M,costheta);
  
  if (funcType == 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 SHFunc::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:	1309
Committed:	Fri Jun 25 20:10:53 2004 UTC (20 years ago) by chrisfen
File size:	1657 byte(s)
Log Message:	Coefficients are now multiplied by the normalization factor - no longer do we need normalization in the visualizer
#	User	Rev	Content
1	gezelter	1291	#include <stdio.h>
2			#include <cmath>
3	gezelter	1290	#include "SHFunc.hpp"
4
5			SHFunc::SHFunc() {
6			}
7
8			double SHFunc::getValueAt(double costheta, double phi) {
9
10	chrisfen	1309	double p, phase;
11	gezelter	1291
12			// associated Legendre polynomial
13			p = LegendreP(L,M,costheta);
14	gezelter	1290
15	gezelter	1291	if (funcType == SH_SIN) {
16			phase = sin((double)M * phi);
17			} else {
18			phase = cos((double)M * phi);
19			}
20
21	chrisfen	1309	return coefficientpphase;
22	gezelter	1290
23			}
24	gezelter	1291	//-----------------------------------------------------------------------------//
25			//
26			// double LegendreP (int l, int m, double x);
27			//
28			// Computes the value of the associated Legendre polynomial P_lm (x)
29			// of order l at a given point.
30			//
31			// Input:
32			// l = degree of the polynomial >= 0
33			// m = parameter satisfying 0 <= m <= l,
34			// x = point in which the computation is performed, range -1 <= x <= 1.
35			// Returns:
36			// value of the polynomial in x
37			//
38			//-----------------------------------------------------------------------------//
39
40			double SHFunc::LegendreP (int l, int m, double x) {
41			// check parameters
42			if (m < 0 \|\| m > l \|\| fabs(x) > 1.0) {
43	gezelter	1300	printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x);
44	gezelter	1291	return NAN;
45			}
46
47			double pmm = 1.0;
48			if (m > 0) {
49			double h = sqrt((1.0-x)*(1.0+x)),
50			f = 1.0;
51			for (int i = 1; i <= m; i++) {
52			pmm = -f h;
53			f += 2.0;
54			}
55			}
56			if (l == m)
57			return pmm;
58			else {
59			double pmmp1 = x * (2 * m + 1) * pmm;
60			if (l == (m+1))
61			return pmmp1;
62			else {
63			double pll = 0.0;
64			for (int ll = m+2; ll <= l; ll++) {
65			pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
66			pmm = pmmp1;
67			pmmp1 = pll;
68			}
69			return pll;
70			}
71			}
72			}
73