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 1309 by chrisfen, Fri Jun 25 20:10:53 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 p, phase;
11 +
12 +  // associated Legendre polynomial
13 +  p = LegendreP(L,M,costheta);
14    
15 +  if (funcType == SH_SIN) {
16 +    phase = sin((double)M * phi);
17 +  } else {
18 +    phase = cos((double)M * phi);
19 +  }
20 +  
21 +  return coefficient*p*phase;
22  
23 + }
24 + //-----------------------------------------------------------------------------//
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 +    printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x);
44 +    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 +

Diff Legend

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