# | Line 119 | Line 119 | void SphereHarm::printShapesFileStart(char name[200], | |
---|---|---|
119 | shapes << "end ShapeInfo\n"; | |
120 | } | |
121 | ||
122 | < | void SphereHarm::printToShapesFile(char name[200], int index){ |
122 | > | void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){ |
123 | ofstream shapes(name, ios::app); | |
124 | ||
125 | biggest = 0.0; | |
# | 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++; |
162 | < | if (fabs(sm) > 0.01 * biggest) nfuncs++; |
161 | > | if (fabs(cm) > tolVal * biggest) nfuncs++; |
162 | > | if (fabs(sm) > tolVal * biggest) nfuncs++; |
163 | } | |
164 | } | |
165 | ||
# | 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) |
196 | > | if (fabs(cm) > tolVal * biggest) |
197 | shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n"; | |
198 | < | if (fabs(sm) > 0.01 * biggest) |
198 | > | if (fabs(sm) > tolVal * biggest) |
199 | shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n"; | |
200 | } | |
201 | } | |
# | 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 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |