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, but
  // we use a lookup table in the factorial code:

  // normalization factor:
  f = sqrt( (2*L+1)/(4.0*M_PI) * Fac(L-M) / Fac(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::Fac (int n) {

  static double facn[31] = {
    1.0,
    1.0,
    2.0,
    6.0,
    24.0,
    120.0,
    720.0,
    5040.0,
    40320.0,
    362880.0,
    3628800.0,
    39916800.0,
    479001600.0,
    6227020800.0,
    87178291200.0,
    1.307674368e12,
    2.0922789888e13,
    3.55687428096e14,
    6.402373705728e15,
    1.21645100408832e17,
    2.43290200817664e18,
    5.109094217170944e19,
    1.12400072777760768e21,
    2.585201673888497664e22,
    6.2044840173323943936e23,
    1.5511210043330985984e25,
    4.03291461126605635584e26,
    1.0888869450418352160768e28,
    3.04888344611713860501504e29,
    8.841761993739701954543616e30,
    2.6525285981219105863630848e32
  };
  
  
  static int nmax = 0;
  static double xmin, xmax;
    
  if (n < 0) {
    printf("factorial of negative integer undefined\n");
    return NAN;
  }
  
  if (n <= 30) return facn[n];
  else {
    printf("n is so large that Fac(n) will overflow\n");
    return NAN;
  }  
}
Revision:	1291
Committed:	Wed Jun 23 22:43:54 2004 UTC (20 years ago) by gezelter
File size:	2806 byte(s)
Log Message:	Added evaluation stuff
#	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			// incredibly inefficient way to get the normalization, but
13			// we use a lookup table in the factorial code:
14
15			// normalization factor:
16			f = sqrt( (2L+1)/(4.0M_PI) * Fac(L-M) / Fac(L+M) );
17			// associated Legendre polynomial
18			p = LegendreP(L,M,costheta);
19	gezelter	1290
20	gezelter	1291	if (funcType == SH_SIN) {
21			phase = sin((double)M * phi);
22			} else {
23			phase = cos((double)M * phi);
24			}
25	gezelter	1290
26	gezelter	1291
27			return coefficientfp*phase;
28	gezelter	1290
29			}
30	gezelter	1291	//-----------------------------------------------------------------------------//
31			//
32			// double LegendreP (int l, int m, double x);
33			//
34			// Computes the value of the associated Legendre polynomial P_lm (x)
35			// of order l at a given point.
36			//
37			// Input:
38			// l = degree of the polynomial >= 0
39			// m = parameter satisfying 0 <= m <= l,
40			// x = point in which the computation is performed, range -1 <= x <= 1.
41			// Returns:
42			// value of the polynomial in x
43			//
44			//-----------------------------------------------------------------------------//
45
46			double SHFunc::LegendreP (int l, int m, double x) {
47			// check parameters
48			if (m < 0 \|\| m > l \|\| fabs(x) > 1.0) {
49			printf("LegendreP got a bad argument\n");
50			return NAN;
51			}
52
53			double pmm = 1.0;
54			if (m > 0) {
55			double h = sqrt((1.0-x)*(1.0+x)),
56			f = 1.0;
57			for (int i = 1; i <= m; i++) {
58			pmm = -f h;
59			f += 2.0;
60			}
61			}
62			if (l == m)
63			return pmm;
64			else {
65			double pmmp1 = x * (2 * m + 1) * pmm;
66			if (l == (m+1))
67			return pmmp1;
68			else {
69			double pll = 0.0;
70			for (int ll = m+2; ll <= l; ll++) {
71			pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
72			pmm = pmmp1;
73			pmmp1 = pll;
74			}
75			return pll;
76			}
77			}
78			}
79
80			double SHFunc::Fac (int n) {
81
82			static double facn[31] = {
83			1.0,
84			1.0,
85			2.0,
86			6.0,
87			24.0,
88			120.0,
89			720.0,
90			5040.0,
91			40320.0,
92			362880.0,
93			3628800.0,
94			39916800.0,
95			479001600.0,
96			6227020800.0,
97			87178291200.0,
98			1.307674368e12,
99			2.0922789888e13,
100			3.55687428096e14,
101			6.402373705728e15,
102			1.21645100408832e17,
103			2.43290200817664e18,
104			5.109094217170944e19,
105			1.12400072777760768e21,
106			2.585201673888497664e22,
107			6.2044840173323943936e23,
108			1.5511210043330985984e25,
109			4.03291461126605635584e26,
110			1.0888869450418352160768e28,
111			3.04888344611713860501504e29,
112			8.841761993739701954543616e30,
113			2.6525285981219105863630848e32
114			};
115
116
117			static int nmax = 0;
118			static double xmin, xmax;
119
120			if (n < 0) {
121			printf("factorial of negative integer undefined\n");
122			return NAN;
123			}
124
125			if (n <= 30) return facn[n];
126			else {
127			printf("n is so large that Fac(n) will overflow\n");
128			return NAN;
129			}
130			}