trunk/SHAPES/SHFunc.cpp

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

SHFunc::SHFunc() {
}

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

  double f, p, phase;

  // incredibly inefficient way to get the normalization

  // normalization factor:
  f = sqrt( (2*L+1)/(4.0*M_PI) * Fact(L-M) / Fact(L+M) );
  // 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*f*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\n");
    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;
    }
  }
}

double SHFunc::Fact(double n) {

  if (n < 0.0) return NAN;
  else {
    if (n < 2.0) return 1.0;
    else
      return n*Fact(n-1.0);
  }
  
}
Revision:	1295
Committed:	Thu Jun 24 15:31:52 2004 UTC (21 years, 4 months ago) by gezelter
File size:	1915 byte(s)
Log Message:	renamed forcer, modified factorial, fixed makefile
#	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	gezelter	1291	double f, p, phase;
11
12	gezelter	1295	// incredibly inefficient way to get the normalization
13	gezelter	1291
14			// normalization factor:
15	gezelter	1295	f = sqrt( (2L+1)/(4.0M_PI) * Fact(L-M) / Fact(L+M) );
16	gezelter	1291	// associated Legendre polynomial
17			p = LegendreP(L,M,costheta);
18	gezelter	1290
19	gezelter	1291	if (funcType == SH_SIN) {
20			phase = sin((double)M * phi);
21			} else {
22			phase = cos((double)M * phi);
23			}
24
25			return coefficientfp*phase;
26	gezelter	1290
27			}
28	gezelter	1291	//-----------------------------------------------------------------------------//
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
44			double SHFunc::LegendreP (int l, int m, double x) {
45			// check parameters
46			if (m < 0 \|\| m > l \|\| fabs(x) > 1.0) {
47			printf("LegendreP got a bad argument\n");
48			return NAN;
49			}
50
51			double pmm = 1.0;
52			if (m > 0) {
53			double h = sqrt((1.0-x)*(1.0+x)),
54			f = 1.0;
55			for (int i = 1; i <= m; i++) {
56			pmm = -f h;
57			f += 2.0;
58			}
59			}
60			if (l == m)
61			return pmm;
62			else {
63			double pmmp1 = x * (2 * m + 1) * pmm;
64			if (l == (m+1))
65			return pmmp1;
66			else {
67			double pll = 0.0;
68			for (int ll = m+2; ll <= l; ll++) {
69			pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
70			pmm = pmmp1;
71			pmmp1 = pll;
72			}
73			return pll;
74			}
75			}
76			}
77
78	gezelter	1295	double SHFunc::Fact(double n) {
79	gezelter	1291
80	gezelter	1295	if (n < 0.0) return NAN;
81			else {
82			if (n < 2.0) return 1.0;
83			else
84			return n*Fact(n-1.0);
85	gezelter	1291	}
86
87			}