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

Comparing trunk/SHAPES/SphereHarm.cpp (file contents):
Revision 1298 by chrisfen, Thu Jun 24 15:47:52 2004 UTC vs.
Revision 1309 by chrisfen, Fri Jun 25 20:10:53 2004 UTC

# Line 130 | Line 130 | void SphereHarm::printToShapesFile(char name[200], int
130        dummy2 = seanindex(-m, l, bw);
131          
132        if (m == 0) {
133 <        cm = rcoeffs[dummy1];
134 <        sm = icoeffs[dummy1];
133 >        cm = normFactor(l,m)*rcoeffs[dummy1];
134 >        sm = normFactor(l,m)*icoeffs[dummy1];
135        } else {
136 <        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
137 <        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
136 >        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1]
137 >                              + rcoeffs[dummy2]);
138 >        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1]
139 >                              - icoeffs[dummy2]);
140        }
141          
142        if (fabs(cm) > biggest) biggest = fabs(cm);
# Line 147 | Line 149 | void SphereHarm::printToShapesFile(char name[200], int
149        dummy2 = seanindex(-m, l, bw);
150          
151        if (m == 0) {
152 <        cm = rcoeffs[dummy1];
153 <        sm = icoeffs[dummy1];
152 >        cm = normFactor(l,m)*rcoeffs[dummy1];
153 >        sm = normFactor(l,m)*icoeffs[dummy1];
154        } else {
155 <        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
156 <        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
155 >        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1]
156 >                              + rcoeffs[dummy2]);
157 >        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1]
158 >                              - icoeffs[dummy2]);
159        }
160          
161        if (fabs(cm) > 0.01 * biggest) nfuncs++;
# Line 180 | Line 184 | void SphereHarm::printToShapesFile(char name[200], int
184        dummy2 = seanindex(-m, l, bw);
185          
186        if (m == 0) {
187 <        cm = rcoeffs[dummy1];
188 <        sm = icoeffs[dummy1];
187 >        cm = normFactor(l,m)*rcoeffs[dummy1];
188 >        sm = normFactor(l,m)*icoeffs[dummy1];
189        } else {
190 <        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
191 <        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
190 >        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1]
191 >                              + rcoeffs[dummy2]);
192 >        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1]
193 >                              - icoeffs[dummy2]);
194        }
195          
196        if (fabs(cm) > 0.01 * biggest)
# Line 206 | Line 212 | void SphereHarm::printToShapesFile(char name[200], int
212    }
213   }
214  
215 + double SphereHarm::normFactor(int L, int M){
216 +  // normalization factor:
217 +  if (L+M > 170){
218 +    printf("Warning: A coefficient was omitted because l + m > 170.\n"
219 +           "\tThe double buffer overflows with factorial calculations\n"
220 +           "\tof 170 and higher.  You should consider using a smaller\n"
221 +           "\tbandwidth if you aren't okay with the loss of the %i, %i\n"
222 +           "\tspherical harmonic.\n", L, M);
223 +    return 0.0;
224 +  }
225 +  else
226 +    return sqrt( (2*L+1)/(4.0*M_PI)*factorialFunc((double)(L-M))
227 +                 / factorialFunc(double(L+M)) );
228 + }
229 +
230 + double SphereHarm::factorialFunc(double n) {
231 +  if (n < 0.0) return NAN;
232 +  else {
233 +    if (n < 2.0) return 1.0;
234 +    else
235 +      return n*factorialFunc(n-1.0);
236 +  }
237 + }

Diff Legend

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