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

vigra/bessel.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 2010-2011 by Ullrich Koethe                  */
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_BESSEL_HXX
00037 #define VIGRA_BESSEL_HXX
00038 
00039 #include "mathutil.hxx"
00040 
00041 #ifdef HasBoostMath
00042 #include <boost/math/special_functions/bessel.hpp>
00043 #endif
00044 
00045 namespace vigra {
00046 
00047 /** \addtogroup MathFunctions 
00048 */
00049 //@{
00050 
00051 namespace detail {
00052 
00053 template <class REAL>
00054 int msta1(REAL x, int mp)
00055 {
00056     double a0,f0,f1,f;
00057     int i,n0,n1,nn;
00058 
00059     a0 = abs(x);
00060     n0 = (int)(1.1*a0)+1;
00061     f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-mp;
00062     n1 = n0+5;
00063     f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-mp;
00064     for(i=0;i<20;i++) 
00065     {
00066         nn = int(n1-(n1-n0)/(1.0-f0/f1));
00067         f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-mp;
00068         if(abs(nn-n1) < 1) 
00069             break;
00070         n0 = n1;
00071         f0 = f1;
00072         n1 = nn;
00073         f1 = f;
00074     }
00075     return nn;
00076 }
00077 
00078 template <class REAL>
00079 int msta2(REAL x, int n, int mp)
00080 {
00081     double a0,ejn,hmp,f0,f1,f,obj;
00082     int i,n0,n1,nn;
00083 
00084     a0 = abs(x);
00085     hmp = 0.5*mp;
00086     ejn = 0.5*std::log10(6.28*n) - n*std::log10(1.36*a0/n);
00087     if (ejn <= hmp) 
00088     {
00089         obj = mp;
00090         n0 = (int)(1.1*a0);
00091         if (n0 < 1) 
00092             n0 = 1;
00093     }
00094     else 
00095     {
00096         obj = hmp+ejn;
00097         n0 = n;
00098     }
00099     f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-obj;
00100     n1 = n0+5;
00101     f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-obj;
00102     for (i=0;i<20;i++) 
00103     {
00104         nn = int(n1-(n1-n0)/(1.0-f0/f1));
00105         f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-obj;
00106         if (abs(nn-n1) < 1) 
00107             break;
00108         n0 = n1;
00109         f0 = f1;
00110         n1 = nn;
00111         f1 = f;
00112     }
00113     return nn+10;
00114 }
00115 
00116 //
00117 //  INPUT:
00118 //  double x    -- argument of Bessel function of 1st and 2nd kind.
00119 //  int n       -- order
00120 //
00121 //  OUPUT:
00122 //
00123 //  int nm      -- highest order actually computed (nm <= n)
00124 //  double jn[] -- Bessel function of 1st kind, orders from 0 to nm
00125 //  double yn[] -- Bessel function of 2nd kind, orders from 0 to nm
00126 //
00127 //  Computes Bessel functions of all order up to 'n' using recurrence
00128 //  relations. If 'nm' < 'n' only 'nm' orders are returned.
00129 //
00130 // code has been adapted from C.R. Bond's implementation
00131 // see http://www.crbond.com/math.htm
00132 //
00133 template <class REAL>
00134 void bessjyn(int n, REAL x,int &nm, double *jn, double *yn)
00135 {
00136     double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
00137     double ec,bs,byk,p0,p1,q0,q1;
00138     static double a[] = {
00139         -0.7031250000000000e-1,
00140          0.1121520996093750,
00141         -0.5725014209747314,
00142          6.074042001273483};
00143     static double b[] = {
00144          0.7324218750000000e-1,
00145         -0.2271080017089844,
00146          1.727727502584457,
00147         -2.438052969955606e1};
00148     static double a1[] = {
00149          0.1171875,
00150         -0.1441955566406250,
00151          0.6765925884246826,
00152         -6.883914268109947};
00153     static double b1[] = {
00154        -0.1025390625,
00155         0.2775764465332031,
00156        -1.993531733751297,
00157         2.724882731126854e1};
00158         
00159     int i,k,m;
00160     nm = n;
00161     if (x < 1e-15) 
00162     {
00163         for (i=0;i<=n;i++) 
00164         {
00165             jn[i] = 0.0;
00166             yn[i] = -1e308;
00167         }
00168         jn[0] = 1.0;
00169         return;
00170     }
00171     if (x <= 300.0 || n > (int)(0.9*x)) 
00172     {
00173         if (n == 0) 
00174             nm = 1;
00175         m = msta1(x,200);
00176         if (m < nm) 
00177             nm = m;
00178         else 
00179             m = msta2(x,nm,15);
00180         bs = 0.0;
00181         su = 0.0;
00182         sv = 0.0;
00183         f2 = 0.0;
00184         f1 = 1.0e-100;
00185         for (k = m;k>=0;k--) 
00186         {
00187             f = 2.0*(k+1.0)/x*f1 - f2;
00188             if (k <= nm) 
00189                 jn[k] = f;
00190             if ((k == 2*(int)(k/2)) && (k != 0)) 
00191             {
00192                 bs += 2.0*f;
00193                 su += (-1)*((k & 2)-1)*f/(double)k;
00194             }
00195             else if (k > 1) 
00196             {
00197                 sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0);
00198             }
00199             f2 = f1;
00200             f1 = f;
00201         }
00202         s0 = bs+f;
00203         for (k=0;k<=nm;k++) 
00204         {
00205             jn[k] /= s0;
00206         }
00207         ec = std::log(0.5*x) + M_EULER_GAMMA;
00208         by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
00209         yn[0] = by0;
00210         by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
00211         yn[1] = by1;
00212     }
00213     else 
00214     {
00215         t1 = x-M_PI_4;
00216         p0 = 1.0;
00217         q0 = -0.125/x;
00218         for (k=0;k<4;k++) 
00219         {
00220             p0 += a[k]*std::pow(x,-2*k-2);
00221             q0 += b[k]*std::pow(x,-2*k-3);
00222         }
00223         cu = std::sqrt(M_2_PI/x);
00224         bj0 = cu*(p0*std::cos(t1)-q0*std::sin(t1));
00225         by0 = cu*(p0*std::sin(t1)+q0*std::cos(t1));
00226         jn[0] = bj0;
00227         yn[0] = by0;
00228         t2 = x-0.75*M_PI;
00229         p1 = 1.0;
00230         q1 = 0.375/x;
00231         for (k=0;k<4;k++) 
00232         {
00233             p1 += a1[k]*std::pow(x,-2*k-2);
00234             q1 += b1[k]*std::pow(x,-2*k-3);
00235         }
00236         bj1 = cu*(p1*std::cos(t2)-q1*std::sin(t2));
00237         by1 = cu*(p1*std::sin(t2)+q1*std::cos(t2));
00238         jn[1] = bj1;
00239         yn[1] = by1;
00240         for (k=2;k<=nm;k++) 
00241         {
00242             bjk = 2.0*(k-1.0)*bj1/x-bj0;
00243             jn[k] = bjk;
00244             bj0 = bj1;
00245             bj1 = bjk;
00246         }
00247     }
00248     for (k=2;k<=nm;k++) 
00249     {
00250         byk = 2.0*(k-1.0)*by1/x-by0;
00251         yn[k] = byk;
00252         by0 = by1;
00253         by1 = byk;
00254     }
00255 }
00256 
00257 
00258  
00259 } // namespace detail
00260 
00261     /*! Bessel function of the first kind. 
00262 
00263         Computes the value of BesselJ of integer order <tt>n</tt> and argument <tt>x</tt>.
00264         Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>.
00265 
00266         This function wraps a number of existing implementations and falls back to 
00267         a rather slow algorithm if none of them is available. In particular,
00268         it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native 
00269         implementations on gcc and MSVC otherwise.
00270 
00271         <b>\#include</b> <vigra/bessel.hxx><br>
00272         Namespace: vigra
00273     */
00274 inline double besselJ(int n, double x)
00275 {
00276     if(x < 0.0)
00277         throw std::domain_error("besselJ(n, x): x cannot be negative");
00278     if(x < 1e-15)
00279         return n == 0 ? 1.0 : 0.0;
00280 #if defined(HasBoostMath)
00281     return boost::math::cyl_bessel_j((double)n, x);
00282 #elif defined(__GNUC__)
00283     return ::jn(n, x);
00284 #elif defined(_MSC_VER)
00285     return _jn(n, x);
00286 #else
00287     int an = abs(n), nr = n, s = an+2;
00288     ArrayVector<double> t(2*s);
00289     detail::bessjyn(an, x, nr, &t[0], &t[s]);
00290     if(n < 0 && odd(an))
00291         return -t[an];
00292     else
00293         return  t[an];
00294 #endif
00295 }
00296 
00297     /*! Bessel function of the second kind. 
00298 
00299         Computes the value of BesselY of integer order <tt>n</tt> and argument <tt>x</tt>.
00300         Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>.
00301 
00302         This function wraps a number of existing implementations and falls back to 
00303         a rather slow algorithm if none of them is available. In particular,
00304         it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native 
00305         implementations on gcc and MSVC otherwise.
00306 
00307         <b>\#include</b> <vigra/bessel.hxx><br>
00308         Namespace: vigra
00309     */
00310 inline double besselY(int n, double x)
00311 {
00312     if(x < 0.0)
00313         throw std::domain_error("besselY(n, x): x cannot be negative");
00314     if(x == 0.0 )
00315         return -std::numeric_limits<double>::infinity();
00316 #if defined(HasBoostMath)
00317     return boost::math::cyl_neumann((double)n, x);
00318 #elif defined(__GNUC__)
00319     return ::yn(n, x);
00320 #elif defined(_MSC_VER)
00321     return _yn(n, x);
00322 #else
00323     int an = abs(n), nr = n, s = an+2;
00324     ArrayVector<double> t(2*s);
00325     detail::bessjyn(an, x, nr, &t[0], &t[s]);
00326     if(n < 0.0 && odd(n))
00327         return -t[an+s];
00328     else
00329         return  t[an+s];
00330 #endif
00331 }
00332 
00333 //@}
00334 
00335 } // namespace vigra
00336 
00337 #endif // VIGRA_BESSEL_HXX

© 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.8.0 (20 Sep 2011)