bessel_IJK_jd

PURPOSE ^

Header source

SYNOPSIS ^

Header source

DESCRIPTION ^

Header source

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 /* Copied from Numerical Recipees modified for use in MEX files*/
0002 
0003 #include <math.h>
0004 #define ACC 40.0 /*Make larger to increase accuracy. */
0005 #define BIGNO 1.0e10
0006 #define BIGNI 1.0e-10
0007 
0008 double bessj0(double x) /* Returns the Bessel function J0 (x) for any real x. */
0009 {
0010     double ax,z;
0011     double xx,y,ans,ans1,ans2; /* Accumulate polynomials in double precision .*/
0012     if ((ax=fabs(x)) < 8.0) { /* Direct rational function fit .*/
0013         y=x*x;
0014         ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
0015             +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
0016         ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
0017             +y*(59272.64853+y*(267.8532712+y*1.0))));
0018         ans=ans1/ans2;
0019     } else { /* Fitting function (6.5.9). */
0020         z=8.0/ax;
0021         y=z*z;
0022         xx=ax-0.785398164;
0023         ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
0024             +y*(-0.2073370639e-5+y*0.2093887211e-6)));
0025         ans2 = -0.1562499995e-1+y*(0.1430488765e-3
0026             +y*(-0.6911147651e-5+y*(0.7621095161e-6
0027             -y*0.934945152e-7)));
0028         ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
0029     }
0030     return ans;
0031 }
0032 
0033 double bessj1(double x) /* Returns the Bessel function J1 (x) for any real x. */
0034 {
0035     double ax,z;
0036     double xx,y,ans,ans1,ans2; /* Accumulate polynomials in double precision. */
0037     if ((ax=fabs(x)) < 8.0) { /* Direct rational approximation. */
0038         y=x*x;
0039         ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
0040             +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
0041         ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
0042             +y*(99447.43394+y*(376.9991397+y*1.0))));
0043         ans=ans1/ans2;
0044     } else { /* Fitting function (6.5.9). */
0045         z=8.0/ax;
0046         y=z*z;
0047         xx=ax-2.356194491;
0048         ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
0049             +y*(0.2457520174e-5+y*(-0.240337019e-6))));
0050         ans2=0.04687499995+y*(-0.2002690873e-3
0051             +y*(0.8449199096e-5+y*(-0.88228987e-6
0052             +y*0.105787412e-6)));
0053         ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
0054         if (x < 0.0) ans = -ans;
0055     }
0056     return ans;
0057 }
0058 
0059 double bessj(int n, double x) /* Returns the Bessel function Jn (x) for any real x and n ≥ 2. */
0060 {
0061     /*double bessj0(double x);
0062     double bessj1(double x);
0063     void nrerror(char error_text[]);*/
0064     int j,jsum,m;
0065     double ax,bj,bjm,bjp,sum,tox,ans;
0066     if (n < 2) mexErrMsgTxt("Index n less than 2 in bessj");
0067     ax=fabs(x);
0068     if (ax == 0.0)
0069         return 0.0;
0070     else if (ax > (double) n) { /* Upwards recurrence from J0 and J1 . */
0071         tox=2.0/ax;
0072         bjm=bessj0(ax);
0073         bj=bessj1(ax);
0074         for (j=1;jelse { /*Downwards recurrence from an even m here computed. */
0081         tox=2.0/ax;                               
0082         m=2*((n+(int) sqrt(ACC*n))/2);
0083         jsum=0; /* jsum will alternate between 0 and 1; when it is 1, we accumulate in sum the even terms in (5.5.16). */
0084         bjp=ans=sum=0.0;                          
0085         bj=1.0;                                   
0086         for (j=m;j>0;j--) { /* The downward recurrence. */
0087              bjm=j*tox*bj-bjp;
0088              bjp=bj;
0089              bj=bjm;
0090              if (fabs(bj) > BIGNO) { /*Renormalize to prevent overflows. */
0091                  bj *= BIGNI;
0092                  bjp *= BIGNI;
0093                  ans *= BIGNI;
0094                  sum *= BIGNI;
0095              }
0096              if (jsum) sum += bj; /*Accumulate the sum. */
0097              jsum=!jsum; /* Change 0 to 1 or vice versa. */
0098              if (j == n) ans=bjp; /*Save the unnormalized answer. */
0099         }
0100         sum=2.0*sum-bj; /*Compute (5.5.16) and use it to normalize the answer. */
0101         ans /= sum;                         
0102     }
0103     return x < 0.0 && (n & 1) ? -ans : ans;
0104 }
0105 
0106 double bessjn(int n, double x) /* Returns the Bessel function Jn (x) for any real x and integer n. */
0107 {
0108     if (n < 0) return pow(-1,-n)*bessjn(-n,x);
0109     else if (n == 0) return bessj0(x);
0110     else if (n == 1) return bessj1(x);
0111     else return bessj(n,x);
0112 }
0113 
0114 double bessi0(double x) /* Returns the modified Bessel function I0(x) for any real x. */
0115 {
0116     double  ax,ans;
0117     double y; /* Accumulate polynomials in double precision. */
0118     if ((ax=fabs(x)) < 3.75) { /* Polynomial fit. */
0119         y=x/3.75;
0120         y*=y;
0121         ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
0122     } else {
0123         y=3.75/ax;
0124         ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1 +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2 +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1 +y*0.392377e-2))))))));
0125     }
0126     return ans;
0127 }
0128 
0129 double bessi1(double x) /* Returns the modified Bessel function I1(x) for any real x. */
0130 {
0131     double ax,ans;
0132     double y; /* Accumulate polynomials in double precision. */
0133     if ((ax=fabs(x)) < 3.75) { /* Polynomial fit. */
0134         y=x/3.75;
0135         y*=y;
0136         ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934 +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
0137     } else {
0138         y=3.75/ax;
0139         ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -y*0.420059e-2));
0140         ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2 +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
0141         ans *= (exp(ax)/sqrt(ax));
0142     }
0143     return x < 0.0 ? -ans : ans;
0144 }
0145 
0146 double bessi(int n, double x) /* Returns the modified Bessel function In(x) for any real x and n ? 2. */
0147 {
0148     /*double bessi0(double x);
0149     void nrerror(char error_text[]);*/
0150     int j;
0151     double bi,bim,bip,tox,ans;
0152     if (n < 2) mexErrMsgTxt("Index n less than 2 in bessi");
0153     if (x == 0.0)
0154         return 0.0;
0155     else {
0156         tox=2.0/fabs(x);
0157         bip=ans=0.0;
0158         bi=1.0;
0159         for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) { /* Downward recurrence from even */
0160             bim=bip+j*tox*bi; /* m. */
0161             bip=bi;
0162             bi=bim;
0163             if (fabs(bi) > BIGNO) { /* Renormalize to prevent overflows. */
0164                 ans *= BIGNI;
0165                 bi *= BIGNI;
0166                 bip *= BIGNI;
0167             }
0168             if (j == n) ans=bip;
0169         }
0170         ans *= bessi0(x)/bi; /* Normalize with bessi0. */
0171         return x < 0.0 && (n & 1) ? -ans : ans;
0172     }
0173 }
0174 
0175 double bessin(int n, double x) /* Returns the modified Bessel function In (x) for any real x and integer n. */
0176 {
0177     if (n < 0) return bessin(-n,x);
0178     else if (n == 0) return bessi0(x);
0179     else if (n == 1) return bessi1(x);
0180     else return bessi(n,x);
0181 }
0182 
0183 double bessk0(double x) /* Returns the modified Bessel function K0(x) for positive real x. */
0184 {
0185     /* double bessi0(double x); */
0186     double y,ans; /* Accumulate polynomials in double precision. */
0187     if (x <= 2.0) { /* Polynomial fit. */
0188         y=x*x/4.0;
0189         ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420 +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2 +y*(0.10750e-3+y*0.74e-5))))));
0190     } else {
0191         y=2.0/x;
0192         ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1 +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2 +y*(-0.251540e-2+y*0.53208e-3))))));
0193         
0194         /* mexPrintf("x,y,ans : %g,%g,%g\n",x,y,ans); */
0195         
0196     }
0197          
0198     return ans;
0199 }
0200 
0201 double bessk1(double x) /* Returns the modified Bessel function K1(x) for positive real x. */
0202     {
0203     /* double bessi1(double x); */
0204     double y,ans; /* Accumulate polynomials in double precision. */
0205     if (x <= 2.0) { /* Polynomial fit. */
0206         y=x*x/4.0;
0207         ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144 +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1 +y*(-0.110404e-2+y*(-0.4686e-4)))))));
0208     } else {
0209         y=2.0/x;
0210         ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619 +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +y*(0.325614e-2+y*(-0.68245e-3)))))));
0211     }
0212     return ans;
0213 }
0214 
0215 double bessk(int n, double x) /* Returns the modified Bessel function Kn(x) for positive x and n >= 2. */
0216     {
0217     /* double bessk0(double x);
0218     double bessk1(double x); 
0219     void nrerror(char error_text[]);*/
0220     int j;
0221     double bk,bkm,bkp,tox;
0222     if (n < 2) mexErrMsgTxt("Index n less than 2 in bessk");
0223     tox=2.0/x;
0224     bkm=bessk0(x); /*Upward recurrence for all x... */
0225     bk=bessk1(x);
0226     for (j=1;j/*...and here it is. */
0227         bkp=bkm+j*tox*bk;
0228         bkm=bk;
0229         bk=bkp;
0230     }
0231     
0232     /* mexPrintf("x,bk0,bk1,bk : %g,%g,%g,%g\n",x,bessk0(x),bessk1(x),bk); */
0233     
0234     return bk;
0235 }
0236 
0237 double besskn(int n, double x) /* Returns the modified Bessel function Kn (x) for any real x and integer n. */
0238 {
0239     if (n < 0) return besskn(-n,x);
0240     else if (n == 0) return bessk0(x);
0241     else if (n == 1) return bessk1(x);
0242     else return bessk(n,x);
0243 }
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.