[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/clebsch-gordan.hxx VIGRA

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)