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 1300 by gezelter, Thu Jun 24 19:06:47 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
13 +
14 +  // normalization factor:
15 +  f = sqrt( (2*L+1)/(4.0*M_PI) * Fact(L-M) / Fact(L+M) );
16 +  // associated Legendre polynomial
17 +  p = LegendreP(L,M,costheta);
18    
19 +  if (funcType == SH_SIN) {
20 +    phase = sin((double)M * phi);
21 +  } else {
22 +    phase = cos((double)M * phi);
23 +  }
24 +  
25 +  return coefficient*f*p*phase;
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  
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: l = %d\tm = %d\tx = %lf\n", l, m, x);
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 + double SHFunc::Fact(double n) {
79 +
80 +  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 +  }
86 +  
87 + }

Diff Legend

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