1 |
|
#include "SphereHarm.hpp" |
2 |
– |
#include "cospmls.h" |
3 |
– |
#include "makeweights.h" |
4 |
– |
#include "FST_semi_memo.h" |
2 |
|
|
3 |
|
SphereHarm::SphereHarm( int bandWidth ){ |
4 |
|
bw = bandWidth; |
109 |
|
weights ); |
110 |
|
} |
111 |
|
|
112 |
< |
void SphereHarm::printShapesFileStart(char name[200], string particle, |
112 |
> |
void SphereHarm::printShapesFileStart(char name[200], char particle[80], |
113 |
|
double mass, double momInert[3][3]){ |
114 |
< |
ofstream shapes(name); |
114 |
> |
ofstream myShapes(name); |
115 |
> |
printShapesStreamStart(myShapes, particle, mass, momInert); |
116 |
> |
} |
117 |
> |
|
118 |
> |
void SphereHarm::printShapesStreamStart(ostream& shapes, char particle[80], |
119 |
> |
double mass, double momInert[3][3]){ |
120 |
> |
shapes << "begin ShapeInfo\n"; |
121 |
> |
shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n"; |
122 |
|
shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t" |
123 |
< |
<< momInert[1][1] << "\t" << momInert[2][2] << "\n\n"; |
123 |
> |
<< momInert[1][1] << "\t" << momInert[2][2] << "\n"; |
124 |
> |
shapes << "end ShapeInfo\n"; |
125 |
|
} |
126 |
|
|
127 |
< |
void SphereHarm::printToShapesFile(char name[200], int index){ |
127 |
> |
|
128 |
> |
|
129 |
> |
void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){ |
130 |
|
ofstream shapes(name, ios::app); |
131 |
+ |
|
132 |
+ |
printToShapesStream(shapes, index, tolVal); |
133 |
+ |
} |
134 |
+ |
|
135 |
+ |
void SphereHarm::printToShapesStream(ostream& shapes, int index, double tolVal) { |
136 |
|
|
137 |
|
biggest = 0.0; |
138 |
|
nfuncs = 0; |
142 |
|
dummy2 = seanindex(-m, l, bw); |
143 |
|
|
144 |
|
if (m == 0) { |
145 |
< |
cm = rcoeffs[dummy1]; |
146 |
< |
sm = icoeffs[dummy1]; |
145 |
> |
cm = normFactor(l,m)*rcoeffs[dummy1]; |
146 |
> |
sm = normFactor(l,m)*icoeffs[dummy1]; |
147 |
|
} else { |
148 |
< |
cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2]; |
149 |
< |
sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2]; |
148 |
> |
cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] |
149 |
> |
+ rcoeffs[dummy2]); |
150 |
> |
sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] |
151 |
> |
- icoeffs[dummy2]); |
152 |
|
} |
153 |
|
|
154 |
|
if (fabs(cm) > biggest) biggest = fabs(cm); |
161 |
|
dummy2 = seanindex(-m, l, bw); |
162 |
|
|
163 |
|
if (m == 0) { |
164 |
< |
cm = rcoeffs[dummy1]; |
165 |
< |
sm = icoeffs[dummy1]; |
164 |
> |
cm = normFactor(l,m)*rcoeffs[dummy1]; |
165 |
> |
sm = normFactor(l,m)*icoeffs[dummy1]; |
166 |
|
} else { |
167 |
< |
cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2]; |
168 |
< |
sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2]; |
167 |
> |
cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] |
168 |
> |
+ rcoeffs[dummy2]); |
169 |
> |
sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] |
170 |
> |
- icoeffs[dummy2]); |
171 |
|
} |
172 |
|
|
173 |
< |
if (fabs(cm) > 0.01 * biggest) nfuncs++; |
174 |
< |
if (fabs(sm) > 0.01 * biggest) nfuncs++; |
173 |
> |
if (fabs(cm) > tolVal * biggest) nfuncs++; |
174 |
> |
if (fabs(sm) > tolVal * biggest) nfuncs++; |
175 |
|
} |
176 |
|
} |
177 |
|
|
196 |
|
dummy2 = seanindex(-m, l, bw); |
197 |
|
|
198 |
|
if (m == 0) { |
199 |
< |
cm = rcoeffs[dummy1]; |
200 |
< |
sm = icoeffs[dummy1]; |
199 |
> |
cm = normFactor(l,m)*rcoeffs[dummy1]; |
200 |
> |
sm = normFactor(l,m)*icoeffs[dummy1]; |
201 |
|
} else { |
202 |
< |
cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2]; |
203 |
< |
sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2]; |
202 |
> |
cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] |
203 |
> |
+ rcoeffs[dummy2]); |
204 |
> |
sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] |
205 |
> |
- icoeffs[dummy2]); |
206 |
|
} |
207 |
|
|
208 |
< |
if (fabs(cm) > 0.01 * biggest) |
209 |
< |
shapes << l << "\t" << m << "\tcos\t" << cm << "\n"; |
210 |
< |
if (fabs(sm) > 0.01 * biggest) |
211 |
< |
shapes << l << "\t" << m << "\tsin\t" << sm << "\n"; |
208 |
> |
if (fabs(cm) > tolVal * biggest) |
209 |
> |
shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n"; |
210 |
> |
if (fabs(sm) > tolVal * biggest) |
211 |
> |
shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n"; |
212 |
|
} |
213 |
|
} |
214 |
|
switch(index){ |
215 |
|
case 0:{ |
216 |
< |
shapes << "\nend ContactFunctions\n"; |
216 |
> |
shapes << "end ContactFunctions\n"; |
217 |
|
}; break; |
218 |
|
case 1:{ |
219 |
< |
shapes << "\nend RangeFunctions\n"; |
219 |
> |
shapes << "end RangeFunctions\n"; |
220 |
|
}; break; |
221 |
|
case 2:{ |
222 |
< |
shapes << "\nend StrengthFunctions\n"; |
222 |
> |
shapes << "end StrengthFunctions\n"; |
223 |
|
}; break; |
224 |
|
} |
225 |
|
} |
226 |
|
|
227 |
+ |
double SphereHarm::normFactor(int L, int M){ |
228 |
+ |
// normalization factor: |
229 |
+ |
if (L+M > 170){ |
230 |
+ |
printf("Warning: A coefficient was omitted because l + m > 170.\n" |
231 |
+ |
"\tThe double buffer overflows with factorial calculations\n" |
232 |
+ |
"\tof 170 and higher. You should consider using a smaller\n" |
233 |
+ |
"\tbandwidth if you aren't okay with the loss of the %i, %i\n" |
234 |
+ |
"\tspherical harmonic.\n", L, M); |
235 |
+ |
return 0.0; |
236 |
+ |
} |
237 |
+ |
else |
238 |
+ |
return sqrt( (2*L+1)/(4.0*M_PI)*factorialFunc((double)(L-M)) |
239 |
+ |
/ factorialFunc(double(L+M)) ); |
240 |
+ |
} |
241 |
+ |
|
242 |
+ |
double SphereHarm::factorialFunc(double n) { |
243 |
+ |
if (n < 0.0) return NAN; |
244 |
+ |
else { |
245 |
+ |
if (n < 2.0) return 1.0; |
246 |
+ |
else |
247 |
+ |
return n*factorialFunc(n-1.0); |
248 |
+ |
} |
249 |
+ |
} |