[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/clebsch-gordan.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2009-2010 by Ullrich Koethe and Janis Fehr */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_CLEBSCH_GORDAN_HXX 00037 #define VIGRA_CLEBSCH_GORDAN_HXX 00038 00039 #include "config.hxx" 00040 #include "numerictraits.hxx" 00041 #include "error.hxx" 00042 #include "mathutil.hxx" 00043 #include "array_vector.hxx" 00044 00045 namespace vigra { 00046 00047 namespace { 00048 00049 void ThreeJSymbolM(double l1, double l2, double l3, double m1, 00050 double &m2min, double &m2max, double *thrcof, int ndim, 00051 int &errflag) 00052 { 00053 ContractViolation err; 00054 const double zero = 0.0, eps = 0.01, one = 1.0, two = 2.0; 00055 00056 int nfin, nlim, i, n, index, lstep, nfinp1, nfinp2, nfinp3, nstep2; 00057 double oldfac, dv, newfac, sumbac = 0.0, thresh, a1s, sumfor, sumuni, 00058 sum1, sum2, x, y, m2, m3, x1, x2, x3, y1, y2, y3, cnorm, 00059 ratio, a1, c1, c2, c1old = 0.0, sign1, sign2; 00060 00061 // Parameter adjustments 00062 --thrcof; 00063 00064 errflag = 0; 00065 00066 // "hugedouble" is the square root of one twentieth of the largest floating 00067 // point number, approximately. 00068 double hugedouble = std::sqrt(NumericTraits<double>::max() / 20.0), 00069 srhuge = std::sqrt(hugedouble), 00070 tiny = one / hugedouble, 00071 srtiny = one / srhuge; 00072 00073 // lmatch = zero 00074 00075 // Check error conditions 1, 2, and 3. 00076 if (l1 - abs(m1) + eps < zero 00077 || std::fmod(l1 + abs(m1) + eps, one) >= eps + eps) 00078 { 00079 errflag = 1; 00080 err << "ThreeJSymbolM: l1-abs(m1) less than zero or l1+abs(m1) not integer.\n"; 00081 throw err; 00082 } 00083 else if (l1+l2-l3 < -eps || l1-l2+l3 < -eps || -(l1) + l2+l3 < -eps) 00084 { 00085 errflag = 2; 00086 err << " ThreeJSymbolM: l1, l2, l3 do not satisfy triangular condition:" 00087 << l1 << " " << l2 << " " << l3 << "\n"; 00088 throw err; 00089 } 00090 else if (std::fmod(l1 + l2 + l3 + eps, one) >= eps + eps) 00091 { 00092 errflag = 3; 00093 err << " ThreeJSymbolM: l1+l2+l3 not integer.\n"; 00094 throw err; 00095 } 00096 00097 // limits for m2 00098 m2min = std::max(-l2,-l3-m1); 00099 m2max = std::min(l2,l3-m1); 00100 00101 // Check error condition 4. 00102 if (std::fmod(m2max - m2min + eps, one) >= eps + eps) { 00103 errflag = 4; 00104 err << " ThreeJSymbolM: m2max-m2min not integer.\n"; 00105 throw err; 00106 } 00107 if (m2min < m2max - eps) 00108 goto L20; 00109 if (m2min < m2max + eps) 00110 goto L10; 00111 00112 // Check error condition 5. 00113 errflag = 5; 00114 err << " ThreeJSymbolM: m2min greater than m2max.\n"; 00115 throw err; 00116 00117 // This is reached in case that m2 and m3 can take only one value. 00118 L10: 00119 // mscale = 0 00120 thrcof[1] = (odd(int(abs(l2-l3-m1)+eps)) 00121 ? -one 00122 : one) / std::sqrt(l1+l2+l3+one); 00123 return; 00124 00125 // This is reached in case that M1 and M2 take more than one value. 00126 L20: 00127 // mscale = 0 00128 nfin = int(m2max - m2min + one + eps); 00129 if (ndim - nfin >= 0) 00130 goto L23; 00131 00132 // Check error condition 6. 00133 00134 errflag = 6; 00135 err << " ThreeJSymbolM: Dimension of result array for 3j coefficients too small.\n"; 00136 throw err; 00137 00138 // Start of forward recursion from m2 = m2min 00139 00140 L23: 00141 m2 = m2min; 00142 thrcof[1] = srtiny; 00143 newfac = 0.0; 00144 c1 = 0.0; 00145 sum1 = tiny; 00146 00147 lstep = 1; 00148 L30: 00149 ++lstep; 00150 m2 += one; 00151 m3 = -m1 - m2; 00152 00153 oldfac = newfac; 00154 a1 = (l2 - m2 + one) * (l2 + m2) * (l3 + m3 + one) * (l3 - m3); 00155 newfac = std::sqrt(a1); 00156 00157 dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one) 00158 - (l2+m2-one) * (l3-m3-one); 00159 00160 if (lstep - 2 > 0) 00161 c1old = abs(c1); 00162 00163 // L32: 00164 c1 = -dv / newfac; 00165 00166 if (lstep > 2) 00167 goto L60; 00168 00169 // If m2 = m2min + 1, the third term in the recursion equation vanishes, 00170 // hence 00171 00172 x = srtiny * c1; 00173 thrcof[2] = x; 00174 sum1 += tiny * c1 * c1; 00175 if (lstep == nfin) 00176 goto L220; 00177 goto L30; 00178 00179 L60: 00180 c2 = -oldfac / newfac; 00181 00182 // Recursion to the next 3j coefficient 00183 x = c1 * thrcof[lstep-1] + c2 * thrcof[lstep-2]; 00184 thrcof[lstep] = x; 00185 sumfor = sum1; 00186 sum1 += x * x; 00187 if (lstep == nfin) 00188 goto L100; 00189 00190 // See if last unnormalized 3j coefficient exceeds srhuge 00191 if (abs(x) < srhuge) 00192 goto L80; 00193 00194 // This is reached if last 3j coefficient larger than srhuge, 00195 // so that the recursion series thrcof(1), ... , thrcof(lstep) 00196 // has to be rescaled to prevent overflow 00197 00198 // mscale = mscale + 1 00199 for (i = 1; i <= lstep; ++i) 00200 { 00201 if (abs(thrcof[i]) < srtiny) 00202 thrcof[i] = zero; 00203 thrcof[i] /= srhuge; 00204 } 00205 sum1 /= hugedouble; 00206 sumfor /= hugedouble; 00207 x /= srhuge; 00208 00209 // As long as abs(c1) is decreasing, the recursion proceeds towards 00210 // increasing 3j values and, hence, is numerically stable. Once 00211 // an increase of abs(c1) is detected, the recursion direction is 00212 // reversed. 00213 00214 L80: 00215 if (c1old - abs(c1) > 0.0) 00216 goto L30; 00217 00218 // Keep three 3j coefficients aroundi mmatch for comparison later 00219 // with backward recursion values. 00220 00221 L100: 00222 // mmatch = m2 - 1 00223 nstep2 = nfin - lstep + 3; 00224 x1 = x; 00225 x2 = thrcof[lstep-1]; 00226 x3 = thrcof[lstep-2]; 00227 00228 // Starting backward recursion from m2max taking nstep2 steps, so 00229 // that forwards and backwards recursion overlap at the three points 00230 // m2 = mmatch+1, mmatch, mmatch-1. 00231 00232 nfinp1 = nfin + 1; 00233 nfinp2 = nfin + 2; 00234 nfinp3 = nfin + 3; 00235 thrcof[nfin] = srtiny; 00236 sum2 = tiny; 00237 00238 m2 = m2max + two; 00239 lstep = 1; 00240 L110: 00241 ++lstep; 00242 m2 -= one; 00243 m3 = -m1 - m2; 00244 oldfac = newfac; 00245 a1s = (l2-m2+two) * (l2+m2-one) * (l3+m3+two) * (l3-m3-one); 00246 newfac = std::sqrt(a1s); 00247 dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one) 00248 - (l2+m2-one) * (l3-m3-one); 00249 c1 = -dv / newfac; 00250 if (lstep > 2) 00251 goto L120; 00252 00253 // if m2 = m2max + 1 the third term in the recursion equation vanishes 00254 00255 y = srtiny * c1; 00256 thrcof[nfin - 1] = y; 00257 if (lstep == nstep2) 00258 goto L200; 00259 sumbac = sum2; 00260 sum2 += y * y; 00261 goto L110; 00262 00263 L120: 00264 c2 = -oldfac / newfac; 00265 00266 // Recursion to the next 3j coefficient 00267 00268 y = c1 * thrcof[nfinp2 - lstep] + c2 * thrcof[nfinp3 - lstep]; 00269 00270 if (lstep == nstep2) 00271 goto L200; 00272 00273 thrcof[nfinp1 - lstep] = y; 00274 sumbac = sum2; 00275 sum2 += y * y; 00276 00277 // See if last 3j coefficient exceeds SRHUGE 00278 00279 if (abs(y) < srhuge) 00280 goto L110; 00281 00282 // This is reached if last 3j coefficient larger than srhuge, 00283 // so that the recursion series thrcof(nfin), ... , thrcof(nfin-lstep+1) 00284 // has to be rescaled to prevent overflow. 00285 00286 // mscale = mscale + 1 00287 for (i = 1; i <= lstep; ++i) 00288 { 00289 index = nfin - i + 1; 00290 if (abs(thrcof[index]) < srtiny) 00291 thrcof[index] = zero; 00292 thrcof[index] /= srhuge; 00293 } 00294 sum2 /= hugedouble; 00295 sumbac /= hugedouble; 00296 00297 goto L110; 00298 00299 // The forward recursion 3j coefficients x1, x2, x3 are to be matched 00300 // with the corresponding backward recursion values y1, y2, y3. 00301 00302 L200: 00303 y3 = y; 00304 y2 = thrcof[nfinp2-lstep]; 00305 y1 = thrcof[nfinp3-lstep]; 00306 00307 // Determine now ratio such that yi = ratio * xi (i=1,2,3) holds 00308 // with minimal error. 00309 00310 ratio = (x1*y1 + x2*y2 + x3*y3) / (x1*x1 + x2*x2 + x3*x3); 00311 nlim = nfin - nstep2 + 1; 00312 00313 if (abs(ratio) < one) 00314 goto L211; 00315 for (n = 1; n <= nlim; ++n) 00316 thrcof[n] = ratio * thrcof[n]; 00317 sumuni = ratio * ratio * sumfor + sumbac; 00318 goto L230; 00319 00320 L211: 00321 ++nlim; 00322 ratio = one / ratio; 00323 for (n = nlim; n <= nfin; ++n) 00324 thrcof[n] = ratio * thrcof[n]; 00325 sumuni = sumfor + ratio * ratio * sumbac; 00326 goto L230; 00327 L220: 00328 sumuni = sum1; 00329 00330 // Normalize 3j coefficients 00331 00332 L230: 00333 cnorm = one / std::sqrt((l1+l1+one) * sumuni); 00334 00335 // Sign convention for last 3j coefficient determines overall phase 00336 00337 sign1 = sign(thrcof[nfin]); 00338 sign2 = odd(int(abs(l2-l3-m1)+eps)) 00339 ? -one 00340 : one; 00341 if (sign1 * sign2 <= 0.0) 00342 goto L235; 00343 else 00344 goto L236; 00345 00346 L235: 00347 cnorm = -cnorm; 00348 00349 L236: 00350 if (abs(cnorm) < one) 00351 goto L250; 00352 00353 for (n = 1; n <= nfin; ++n) 00354 thrcof[n] = cnorm * thrcof[n]; 00355 return; 00356 00357 L250: 00358 thresh = tiny / abs(cnorm); 00359 for (n = 1; n <= nfin; ++n) 00360 { 00361 if (abs(thrcof[n]) < thresh) 00362 thrcof[n] = zero; 00363 thrcof[n] = cnorm * thrcof[n]; 00364 } 00365 } 00366 00367 } // anonymous namespace 00368 00369 inline 00370 double clebschGordan (double l1, double m1, double l2, double m2, double l3, double m3) 00371 { 00372 const double err = 0.01; 00373 double CG = 0.0, m2min, m2max, *cofp; 00374 // static array for calculation of 3-j symbols 00375 const int ncof = 100; 00376 double cof[ncof]; 00377 // reset error flag 00378 int errflag = 0; 00379 ContractViolation Err; 00380 00381 // Check for physical restriction. 00382 // All other restrictions are checked by the 3-j symbol routine. 00383 if ( abs(m1 + m2 - m3) > err) 00384 { 00385 errflag = 7; 00386 Err << " clebschGordan: m1 + m2 - m3 is not zero.\n"; 00387 throw Err; 00388 } 00389 // calculate minimum storage size needed for ThreeJSymbolM() 00390 // if the dimension becomes negative the 3-j routine will capture it 00391 int njm = roundi(std::min(l2,l3-m1) - std::max(-l2,-l3-m1) + 1); 00392 00393 // allocate dynamic memory if necessary 00394 ArrayVector<double> cofa; 00395 if(njm > ncof) 00396 { 00397 cofa.resize(njm); 00398 cofp = cofa.begin(); 00399 } 00400 else 00401 { 00402 cofp = cof; 00403 } 00404 00405 // calculate series of 3-j symbols 00406 ThreeJSymbolM (l1,l2,l3,m1, m2min,m2max, cofp,njm, errflag); 00407 // calculated Clebsch-Gordan coefficient 00408 if (! errflag) 00409 { 00410 CG = cofp[roundi(m2-m2min)] * (odd(roundi(l1-l2+m3)) ? -1 : 1) * std::sqrt(2*l3+1); 00411 } 00412 else 00413 { 00414 Err << " clebschGordan: 3jM-sym error.\n"; 00415 throw Err; 00416 } 00417 return CG; 00418 } 00419 00420 } // namespace vigra 00421 00422 #endif // VIGRA_CLEBSCH_GORDAN_HXX 00423
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|