ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/SHAPES/SHFunc.cpp
(Generate patch)

Comparing trunk/SHAPES/SHFunc.cpp (file contents):
Revision 1290 by gezelter, Wed Jun 23 20:53:17 2004 UTC vs.
Revision 1291 by gezelter, Wed Jun 23 22:43:54 2004 UTC

# Line 1 | Line 1
1 + #include <stdio.h>
2 + #include <cmath>
3   #include "SHFunc.hpp"
4  
5   SHFunc::SHFunc() {
# Line 5 | Line 7 | double SHFunc::getValueAt(double costheta, double phi)
7  
8   double SHFunc::getValueAt(double costheta, double phi) {
9  
10 +  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( (2*L+1)/(4.0*M_PI) * Fac(L-M) / Fac(L+M) );
17 +  // associated Legendre polynomial
18 +  p = LegendreP(L,M,costheta);
19    
20 +  if (funcType == SH_SIN) {
21 +    phase = sin((double)M * phi);
22 +  } else {
23 +    phase = cos((double)M * phi);
24 +  }
25  
26 +  
27 +  return coefficient*f*p*phase;
28  
29   }
30 + //-----------------------------------------------------------------------------//
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 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines