vigra - vigra/clebsch-gordan.hxx Source File

[ 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)