1 /************************************************************************/
2 /* */
3 /* Copyright 2009-2010 by Ullrich Koethe and Janis Fehr */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35
36 #ifndef VIGRA_CLEBSCH_GORDAN_HXX
37 #define VIGRA_CLEBSCH_GORDAN_HXX
38
39 #include "config.hxx"
40 #include "numerictraits.hxx"
41 #include "error.hxx"
42 #include "mathutil.hxx"
43 #include "array_vector.hxx"
44
45 namespace vigra {
46
47 namespace {
48
ThreeJSymbolM(double l1,double l2,double l3,double m1,double & m2min,double & m2max,double * thrcof,int ndim,int & errflag)49 void ThreeJSymbolM(double l1, double l2, double l3, double m1,
50 double &m2min, double &m2max, double *thrcof, int ndim,
51 int &errflag)
52 {
53 ContractViolation err;
54 const double zero = 0.0, eps = 0.01, one = 1.0, two = 2.0;
55
56 int nfin, nlim, i, n, index, lstep, nfinp1, nfinp2, nfinp3, nstep2;
57 double oldfac, dv, newfac, sumbac = 0.0, thresh, a1s, sumfor, sumuni,
58 sum1, sum2, x, y, m2, m3, x1, x2, x3, y1, y2, y3, cnorm,
59 ratio, a1, c1, c2, c1old = 0.0, sign1, sign2;
60
61 // Parameter adjustments
62 --thrcof;
63
64 errflag = 0;
65
66 // "hugedouble" is the square root of one twentieth of the largest floating
67 // point number, approximately.
68 double hugedouble = std::sqrt(NumericTraits<double>::max() / 20.0),
69 srhuge = std::sqrt(hugedouble),
70 tiny = one / hugedouble,
71 srtiny = one / srhuge;
72
73 // lmatch = zero
74
75 // Check error conditions 1, 2, and 3.
76 if (l1 - abs(m1) + eps < zero
77 || std::fmod(l1 + abs(m1) + eps, one) >= eps + eps)
78 {
79 errflag = 1;
80 err << "ThreeJSymbolM: l1-abs(m1) less than zero or l1+abs(m1) not integer.\n";
81 throw err;
82 }
83 else if (l1+l2-l3 < -eps || l1-l2+l3 < -eps || -(l1) + l2+l3 < -eps)
84 {
85 errflag = 2;
86 err << " ThreeJSymbolM: l1, l2, l3 do not satisfy triangular condition:"
87 << l1 << " " << l2 << " " << l3 << "\n";
88 throw err;
89 }
90 else if (std::fmod(l1 + l2 + l3 + eps, one) >= eps + eps)
91 {
92 errflag = 3;
93 err << " ThreeJSymbolM: l1+l2+l3 not integer.\n";
94 throw err;
95 }
96
97 // limits for m2
98 m2min = std::max(-l2,-l3-m1);
99 m2max = std::min(l2,l3-m1);
100
101 // Check error condition 4.
102 if (std::fmod(m2max - m2min + eps, one) >= eps + eps) {
103 errflag = 4;
104 err << " ThreeJSymbolM: m2max-m2min not integer.\n";
105 throw err;
106 }
107 if (m2min < m2max - eps)
108 goto L20;
109 if (m2min < m2max + eps)
110 goto L10;
111
112 // Check error condition 5.
113 errflag = 5;
114 err << " ThreeJSymbolM: m2min greater than m2max.\n";
115 throw err;
116
117 // This is reached in case that m2 and m3 can take only one value.
118 L10:
119 // mscale = 0
120 thrcof[1] = (odd(int(abs(l2-l3-m1)+eps))
121 ? -one
122 : one) / std::sqrt(l1+l2+l3+one);
123 return;
124
125 // This is reached in case that M1 and M2 take more than one value.
126 L20:
127 // mscale = 0
128 nfin = int(m2max - m2min + one + eps);
129 if (ndim - nfin >= 0)
130 goto L23;
131
132 // Check error condition 6.
133
134 errflag = 6;
135 err << " ThreeJSymbolM: Dimension of result array for 3j coefficients too small.\n";
136 throw err;
137
138 // Start of forward recursion from m2 = m2min
139
140 L23:
141 m2 = m2min;
142 thrcof[1] = srtiny;
143 newfac = 0.0;
144 c1 = 0.0;
145 sum1 = tiny;
146
147 lstep = 1;
148 L30:
149 ++lstep;
150 m2 += one;
151 m3 = -m1 - m2;
152
153 oldfac = newfac;
154 a1 = (l2 - m2 + one) * (l2 + m2) * (l3 + m3 + one) * (l3 - m3);
155 newfac = std::sqrt(a1);
156
157 dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one)
158 - (l2+m2-one) * (l3-m3-one);
159
160 if (lstep - 2 > 0)
161 c1old = abs(c1);
162
163 // L32:
164 c1 = -dv / newfac;
165
166 if (lstep > 2)
167 goto L60;
168
169 // If m2 = m2min + 1, the third term in the recursion equation vanishes,
170 // hence
171
172 x = srtiny * c1;
173 thrcof[2] = x;
174 sum1 += tiny * c1 * c1;
175 if (lstep == nfin)
176 goto L220;
177 goto L30;
178
179 L60:
180 c2 = -oldfac / newfac;
181
182 // Recursion to the next 3j coefficient
183 x = c1 * thrcof[lstep-1] + c2 * thrcof[lstep-2];
184 thrcof[lstep] = x;
185 sumfor = sum1;
186 sum1 += x * x;
187 if (lstep == nfin)
188 goto L100;
189
190 // See if last unnormalized 3j coefficient exceeds srhuge
191 if (abs(x) < srhuge)
192 goto L80;
193
194 // This is reached if last 3j coefficient larger than srhuge,
195 // so that the recursion series thrcof(1), ... , thrcof(lstep)
196 // has to be rescaled to prevent overflow
197
198 // mscale = mscale + 1
199 for (i = 1; i <= lstep; ++i)
200 {
201 if (abs(thrcof[i]) < srtiny)
202 thrcof[i] = zero;
203 thrcof[i] /= srhuge;
204 }
205 sum1 /= hugedouble;
206 sumfor /= hugedouble;
207 x /= srhuge;
208
209 // As long as abs(c1) is decreasing, the recursion proceeds towards
210 // increasing 3j values and, hence, is numerically stable. Once
211 // an increase of abs(c1) is detected, the recursion direction is
212 // reversed.
213
214 L80:
215 if (c1old - abs(c1) > 0.0)
216 goto L30;
217
218 // Keep three 3j coefficients aroundi mmatch for comparison later
219 // with backward recursion values.
220
221 L100:
222 // mmatch = m2 - 1
223 nstep2 = nfin - lstep + 3;
224 x1 = x;
225 x2 = thrcof[lstep-1];
226 x3 = thrcof[lstep-2];
227
228 // Starting backward recursion from m2max taking nstep2 steps, so
229 // that forwards and backwards recursion overlap at the three points
230 // m2 = mmatch+1, mmatch, mmatch-1.
231
232 nfinp1 = nfin + 1;
233 nfinp2 = nfin + 2;
234 nfinp3 = nfin + 3;
235 thrcof[nfin] = srtiny;
236 sum2 = tiny;
237
238 m2 = m2max + two;
239 lstep = 1;
240 L110:
241 ++lstep;
242 m2 -= one;
243 m3 = -m1 - m2;
244 oldfac = newfac;
245 a1s = (l2-m2+two) * (l2+m2-one) * (l3+m3+two) * (l3-m3-one);
246 newfac = std::sqrt(a1s);
247 dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one)
248 - (l2+m2-one) * (l3-m3-one);
249 c1 = -dv / newfac;
250 if (lstep > 2)
251 goto L120;
252
253 // if m2 = m2max + 1 the third term in the recursion equation vanishes
254
255 y = srtiny * c1;
256 thrcof[nfin - 1] = y;
257 if (lstep == nstep2)
258 goto L200;
259 sumbac = sum2;
260 sum2 += y * y;
261 goto L110;
262
263 L120:
264 c2 = -oldfac / newfac;
265
266 // Recursion to the next 3j coefficient
267
268 y = c1 * thrcof[nfinp2 - lstep] + c2 * thrcof[nfinp3 - lstep];
269
270 if (lstep == nstep2)
271 goto L200;
272
273 thrcof[nfinp1 - lstep] = y;
274 sumbac = sum2;
275 sum2 += y * y;
276
277 // See if last 3j coefficient exceeds SRHUGE
278
279 if (abs(y) < srhuge)
280 goto L110;
281
282 // This is reached if last 3j coefficient larger than srhuge,
283 // so that the recursion series thrcof(nfin), ... , thrcof(nfin-lstep+1)
284 // has to be rescaled to prevent overflow.
285
286 // mscale = mscale + 1
287 for (i = 1; i <= lstep; ++i)
288 {
289 index = nfin - i + 1;
290 if (abs(thrcof[index]) < srtiny)
291 thrcof[index] = zero;
292 thrcof[index] /= srhuge;
293 }
294 sum2 /= hugedouble;
295 sumbac /= hugedouble;
296
297 goto L110;
298
299 // The forward recursion 3j coefficients x1, x2, x3 are to be matched
300 // with the corresponding backward recursion values y1, y2, y3.
301
302 L200:
303 y3 = y;
304 y2 = thrcof[nfinp2-lstep];
305 y1 = thrcof[nfinp3-lstep];
306
307 // Determine now ratio such that yi = ratio * xi (i=1,2,3) holds
308 // with minimal error.
309
310 ratio = (x1*y1 + x2*y2 + x3*y3) / (x1*x1 + x2*x2 + x3*x3);
311 nlim = nfin - nstep2 + 1;
312
313 if (abs(ratio) < one)
314 goto L211;
315 for (n = 1; n <= nlim; ++n)
316 thrcof[n] = ratio * thrcof[n];
317 sumuni = ratio * ratio * sumfor + sumbac;
318 goto L230;
319
320 L211:
321 ++nlim;
322 ratio = one / ratio;
323 for (n = nlim; n <= nfin; ++n)
324 thrcof[n] = ratio * thrcof[n];
325 sumuni = sumfor + ratio * ratio * sumbac;
326 goto L230;
327 L220:
328 sumuni = sum1;
329
330 // Normalize 3j coefficients
331
332 L230:
333 cnorm = one / std::sqrt((l1+l1+one) * sumuni);
334
335 // Sign convention for last 3j coefficient determines overall phase
336
337 sign1 = sign(thrcof[nfin]);
338 sign2 = odd(int(abs(l2-l3-m1)+eps))
339 ? -one
340 : one;
341 if (sign1 * sign2 <= 0.0)
342 goto L235;
343 else
344 goto L236;
345
346 L235:
347 cnorm = -cnorm;
348
349 L236:
350 if (abs(cnorm) < one)
351 goto L250;
352
353 for (n = 1; n <= nfin; ++n)
354 thrcof[n] = cnorm * thrcof[n];
355 return;
356
357 L250:
358 thresh = tiny / abs(cnorm);
359 for (n = 1; n <= nfin; ++n)
360 {
361 if (abs(thrcof[n]) < thresh)
362 thrcof[n] = zero;
363 thrcof[n] = cnorm * thrcof[n];
364 }
365 }
366
367 } // anonymous namespace
368
369 inline
clebschGordan(double l1,double m1,double l2,double m2,double l3,double m3)370 double clebschGordan (double l1, double m1, double l2, double m2, double l3, double m3)
371 {
372 const double err = 0.01;
373 double CG = 0.0, m2min, m2max, *cofp;
374 // array for calculation of 3-j symbols
375 const int ncof = 100;
376 double cof[ncof];
377 // reset error flag
378 int errflag = 0;
379 ContractViolation Err;
380
381 // Check for physical restriction.
382 // All other restrictions are checked by the 3-j symbol routine.
383 if ( abs(m1 + m2 - m3) > err)
384 {
385 errflag = 7;
386 Err << " clebschGordan: m1 + m2 - m3 is not zero.\n";
387 throw Err;
388 }
389 // calculate minimum storage size needed for ThreeJSymbolM()
390 // if the dimension becomes negative the 3-j routine will capture it
391 int njm = roundi(std::min(l2,l3-m1) - std::max(-l2,-l3-m1) + 1);
392
393 // allocate dynamic memory if necessary
394 ArrayVector<double> cofa;
395 if(njm > ncof)
396 {
397 cofa.resize(njm);
398 cofp = cofa.begin();
399 }
400 else
401 {
402 cofp = cof;
403 }
404
405 // calculate series of 3-j symbols
406 ThreeJSymbolM (l1,l2,l3,m1, m2min,m2max, cofp,njm, errflag);
407 // calculated Clebsch-Gordan coefficient
408 if (! errflag)
409 {
410 CG = cofp[roundi(m2-m2min)] * (odd(roundi(l1-l2+m3)) ? -1 : 1) * std::sqrt(2*l3+1);
411 }
412 else
413 {
414 Err << " clebschGordan: 3jM-sym error.\n";
415 throw Err;
416 }
417 return CG;
418 }
419
420 } // namespace vigra
421
422 #endif // VIGRA_CLEBSCH_GORDAN_HXX
423
424