/************************************************************************* Cephes Math Library Release 2.8: June, 2000 Copyright by Stephen L. Moshier Contributors: * Sergey Bochkanov (ALGLIB project). Translation from C to pseudocode. See subroutines comments for additional copyrights. >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; namespace alglib { public class bessel { /************************************************************************* Bessel function of order zero Returns Bessel function of order zero of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used: 2 2 (w - r ) (w - r ) P (w) / Q (w) 1 2 3 8 2 where w = x and the two r's are zeros of the function. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 60000 4.2e-16 1.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselj0(double x) { double result = 0; double xsq = 0; double nn = 0; double pzero = 0; double qzero = 0; double p1 = 0; double q1 = 0; if( (double)(x)<(double)(0) ) { x = -x; } if( (double)(x)>(double)(8.0) ) { besselasympt0(x, ref pzero, ref qzero); nn = x-Math.PI/4; result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn)); return result; } xsq = AP.Math.Sqr(x); p1 = 26857.86856980014981415848441; p1 = -40504123.71833132706360663322+xsq*p1; p1 = 25071582855.36881945555156435+xsq*p1; p1 = -8085222034853.793871199468171+xsq*p1; p1 = 1434354939140344.111664316553+xsq*p1; p1 = -136762035308817138.6865416609+xsq*p1; p1 = 6382059341072356562.289432465+xsq*p1; p1 = -117915762910761053603.8440800+xsq*p1; p1 = 493378725179413356181.6813446+xsq*p1; q1 = 1.0; q1 = 1363.063652328970604442810507+xsq*q1; q1 = 1114636.098462985378182402543+xsq*q1; q1 = 669998767.2982239671814028660+xsq*q1; q1 = 312304311494.1213172572469442+xsq*q1; q1 = 112775673967979.8507056031594+xsq*q1; q1 = 30246356167094626.98627330784+xsq*q1; q1 = 5428918384092285160.200195092+xsq*q1; q1 = 493378725179413356211.3278438+xsq*q1; result = p1/q1; return result; } /************************************************************************* Bessel function of order one Returns Bessel function of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 30000 2.6e-16 1.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselj1(double x) { double result = 0; double s = 0; double xsq = 0; double nn = 0; double pzero = 0; double qzero = 0; double p1 = 0; double q1 = 0; s = Math.Sign(x); if( (double)(x)<(double)(0) ) { x = -x; } if( (double)(x)>(double)(8.0) ) { besselasympt1(x, ref pzero, ref qzero); nn = x-3*Math.PI/4; result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn)); if( (double)(s)<(double)(0) ) { result = -result; } return result; } xsq = AP.Math.Sqr(x); p1 = 2701.122710892323414856790990; p1 = -4695753.530642995859767162166+xsq*p1; p1 = 3413234182.301700539091292655+xsq*p1; p1 = -1322983480332.126453125473247+xsq*p1; p1 = 290879526383477.5409737601689+xsq*p1; p1 = -35888175699101060.50743641413+xsq*p1; p1 = 2316433580634002297.931815435+xsq*p1; p1 = -66721065689249162980.20941484+xsq*p1; p1 = 581199354001606143928.050809+xsq*p1; q1 = 1.0; q1 = 1606.931573481487801970916749+xsq*q1; q1 = 1501793.594998585505921097578+xsq*q1; q1 = 1013863514.358673989967045588+xsq*q1; q1 = 524371026216.7649715406728642+xsq*q1; q1 = 208166122130760.7351240184229+xsq*q1; q1 = 60920613989175217.46105196863+xsq*q1; q1 = 11857707121903209998.37113348+xsq*q1; q1 = 1162398708003212287858.529400+xsq*q1; result = s*x*p1/q1; return result; } /************************************************************************* Bessel function of integer order Returns Bessel function of order n, where n is a (possibly negative) integer. The ratio of jn(x) to j0(x) is computed by backward recurrence. First the ratio jn/jn-1 is found by a continued fraction expansion. Then the recurrence relating successive orders is applied until j0 or j1 is reached. If n = 0 or 1 the routine for j0 or j1 is called directly. ACCURACY: Absolute error: arithmetic range # trials peak rms IEEE 0, 30 5000 4.4e-16 7.9e-17 Not suitable for large n or x. Use jv() (fractional order) instead. Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besseljn(int n, double x) { double result = 0; double pkm2 = 0; double pkm1 = 0; double pk = 0; double xk = 0; double r = 0; double ans = 0; int k = 0; int sg = 0; if( n<0 ) { n = -n; if( n%2==0 ) { sg = 1; } else { sg = -1; } } else { sg = 1; } if( (double)(x)<(double)(0) ) { if( n%2!=0 ) { sg = -sg; } x = -x; } if( n==0 ) { result = sg*besselj0(x); return result; } if( n==1 ) { result = sg*besselj1(x); return result; } if( n==2 ) { if( (double)(x)==(double)(0) ) { result = 0; } else { result = sg*(2.0*besselj1(x)/x-besselj0(x)); } return result; } if( (double)(x)<(double)(AP.Math.MachineEpsilon) ) { result = 0; return result; } k = 53; pk = 2*(n+k); ans = pk; xk = x*x; do { pk = pk-2.0; ans = pk-xk/ans; k = k-1; } while( k!=0 ); ans = x/ans; pk = 1.0; pkm1 = 1.0/ans; k = n-1; r = 2*k; do { pkm2 = (pkm1*r-pk*x)/x; pk = pkm1; pkm1 = pkm2; r = r-2.0; k = k-1; } while( k!=0 ); if( (double)(Math.Abs(pk))>(double)(Math.Abs(pkm1)) ) { ans = besselj1(x)/pk; } else { ans = besselj0(x)/pkm1; } result = sg*ans; return result; } /************************************************************************* Bessel function of the second kind, order zero Returns Bessel function of the second kind, of order zero, of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation R(x) is employed to compute y0(x) = R(x) + 2 * log(x) * j0(x) / PI. Thus a call to j0() is required. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error, when y0(x) < 1; else relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.3e-15 1.6e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ public static double bessely0(double x) { double result = 0; double nn = 0; double xsq = 0; double pzero = 0; double qzero = 0; double p4 = 0; double q4 = 0; if( (double)(x)>(double)(8.0) ) { besselasympt0(x, ref pzero, ref qzero); nn = x-Math.PI/4; result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn)); return result; } xsq = AP.Math.Sqr(x); p4 = -41370.35497933148554125235152; p4 = 59152134.65686889654273830069+xsq*p4; p4 = -34363712229.79040378171030138+xsq*p4; p4 = 10255208596863.94284509167421+xsq*p4; p4 = -1648605817185729.473122082537+xsq*p4; p4 = 137562431639934407.8571335453+xsq*p4; p4 = -5247065581112764941.297350814+xsq*p4; p4 = 65874732757195549259.99402049+xsq*p4; p4 = -27502866786291095837.01933175+xsq*p4; q4 = 1.0; q4 = 1282.452772478993804176329391+xsq*q4; q4 = 1001702.641288906265666651753+xsq*q4; q4 = 579512264.0700729537480087915+xsq*q4; q4 = 261306575504.1081249568482092+xsq*q4; q4 = 91620380340751.85262489147968+xsq*q4; q4 = 23928830434997818.57439356652+xsq*q4; q4 = 4192417043410839973.904769661+xsq*q4; q4 = 372645883898616588198.9980+xsq*q4; result = p4/q4+2/Math.PI*besselj0(x)*Math.Log(x); return result; } /************************************************************************* Bessel function of second kind of order one Returns Bessel function of the second kind of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and a call to j1() is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.0e-15 1.3e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ public static double bessely1(double x) { double result = 0; double nn = 0; double xsq = 0; double pzero = 0; double qzero = 0; double p4 = 0; double q4 = 0; if( (double)(x)>(double)(8.0) ) { besselasympt1(x, ref pzero, ref qzero); nn = x-3*Math.PI/4; result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn)); return result; } xsq = AP.Math.Sqr(x); p4 = -2108847.540133123652824139923; p4 = 3639488548.124002058278999428+xsq*p4; p4 = -2580681702194.450950541426399+xsq*p4; p4 = 956993023992168.3481121552788+xsq*p4; p4 = -196588746272214065.8820322248+xsq*p4; p4 = 21931073399177975921.11427556+xsq*p4; p4 = -1212297555414509577913.561535+xsq*p4; p4 = 26554738314348543268942.48968+xsq*p4; p4 = -99637534243069222259967.44354+xsq*p4; q4 = 1.0; q4 = 1612.361029677000859332072312+xsq*q4; q4 = 1563282.754899580604737366452+xsq*q4; q4 = 1128686837.169442121732366891+xsq*q4; q4 = 646534088126.5275571961681500+xsq*q4; q4 = 297663212564727.6729292742282+xsq*q4; q4 = 108225825940881955.2553850180+xsq*q4; q4 = 29549879358971486742.90758119+xsq*q4; q4 = 5435310377188854170800.653097+xsq*q4; q4 = 508206736694124324531442.4152+xsq*q4; result = x*p4/q4+2/Math.PI*(besselj1(x)*Math.Log(x)-1/x); return result; } /************************************************************************* Bessel function of second kind of integer order Returns Bessel function of order n, where n is a (possibly negative) integer. The function is evaluated by forward recurrence on n, starting with values computed by the routines y0() and y1(). If n = 0 or 1 the routine for y0 or y1 is called directly. ACCURACY: Absolute error, except relative when y > 1: arithmetic domain # trials peak rms IEEE 0, 30 30000 3.4e-15 4.3e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselyn(int n, double x) { double result = 0; int i = 0; double a = 0; double b = 0; double tmp = 0; double s = 0; s = 1; if( n<0 ) { n = -n; if( n%2!=0 ) { s = -1; } } if( n==0 ) { result = bessely0(x); return result; } if( n==1 ) { result = s*bessely1(x); return result; } a = bessely0(x); b = bessely1(x); for(i=1; i<=n-1; i++) { tmp = b; b = 2*i/x*b-a; a = tmp; } result = s*b; return result; } /************************************************************************* Modified Bessel function of order zero Returns modified Bessel function of order zero of the argument. The function is defined as i0(x) = j0( ix ). The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,30 30000 5.8e-16 1.4e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besseli0(double x) { double result = 0; double y = 0; double v = 0; double z = 0; double b0 = 0; double b1 = 0; double b2 = 0; if( (double)(x)<(double)(0) ) { x = -x; } if( (double)(x)<=(double)(8.0) ) { y = x/2.0-2.0; besselmfirstcheb(-4.41534164647933937950E-18, ref b0, ref b1, ref b2); besselmnextcheb(y, 3.33079451882223809783E-17, ref b0, ref b1, ref b2); besselmnextcheb(y, -2.43127984654795469359E-16, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.71539128555513303061E-15, ref b0, ref b1, ref b2); besselmnextcheb(y, -1.16853328779934516808E-14, ref b0, ref b1, ref b2); besselmnextcheb(y, 7.67618549860493561688E-14, ref b0, ref b1, ref b2); besselmnextcheb(y, -4.85644678311192946090E-13, ref b0, ref b1, ref b2); besselmnextcheb(y, 2.95505266312963983461E-12, ref b0, ref b1, ref b2); besselmnextcheb(y, -1.72682629144155570723E-11, ref b0, ref b1, ref b2); besselmnextcheb(y, 9.67580903537323691224E-11, ref b0, ref b1, ref b2); besselmnextcheb(y, -5.18979560163526290666E-10, ref b0, ref b1, ref b2); besselmnextcheb(y, 2.65982372468238665035E-9, ref b0, ref b1, ref b2); besselmnextcheb(y, -1.30002500998624804212E-8, ref b0, ref b1, ref b2); besselmnextcheb(y, 6.04699502254191894932E-8, ref b0, ref b1, ref b2); besselmnextcheb(y, -2.67079385394061173391E-7, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.11738753912010371815E-6, ref b0, ref b1, ref b2); besselmnextcheb(y, -4.41673835845875056359E-6, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.64484480707288970893E-5, ref b0, ref b1, ref b2); besselmnextcheb(y, -5.75419501008210370398E-5, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.88502885095841655729E-4, ref b0, ref b1, ref b2); besselmnextcheb(y, -5.76375574538582365885E-4, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.63947561694133579842E-3, ref b0, ref b1, ref b2); besselmnextcheb(y, -4.32430999505057594430E-3, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.05464603945949983183E-2, ref b0, ref b1, ref b2); besselmnextcheb(y, -2.37374148058994688156E-2, ref b0, ref b1, ref b2); besselmnextcheb(y, 4.93052842396707084878E-2, ref b0, ref b1, ref b2); besselmnextcheb(y, -9.49010970480476444210E-2, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.71620901522208775349E-1, ref b0, ref b1, ref b2); besselmnextcheb(y, -3.04682672343198398683E-1, ref b0, ref b1, ref b2); besselmnextcheb(y, 6.76795274409476084995E-1, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); result = Math.Exp(x)*v; return result; } z = 32.0/x-2.0; besselmfirstcheb(-7.23318048787475395456E-18, ref b0, ref b1, ref b2); besselmnextcheb(z, -4.83050448594418207126E-18, ref b0, ref b1, ref b2); besselmnextcheb(z, 4.46562142029675999901E-17, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.46122286769746109310E-17, ref b0, ref b1, ref b2); besselmnextcheb(z, -2.82762398051658348494E-16, ref b0, ref b1, ref b2); besselmnextcheb(z, -3.42548561967721913462E-16, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.77256013305652638360E-15, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.81168066935262242075E-15, ref b0, ref b1, ref b2); besselmnextcheb(z, -9.55484669882830764870E-15, ref b0, ref b1, ref b2); besselmnextcheb(z, -4.15056934728722208663E-14, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.54008621752140982691E-14, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.85277838274214270114E-13, ref b0, ref b1, ref b2); besselmnextcheb(z, 7.18012445138366623367E-13, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.79417853150680611778E-12, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.32158118404477131188E-11, ref b0, ref b1, ref b2); besselmnextcheb(z, -3.14991652796324136454E-11, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.18891471078464383424E-11, ref b0, ref b1, ref b2); besselmnextcheb(z, 4.94060238822496958910E-10, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.39623202570838634515E-9, ref b0, ref b1, ref b2); besselmnextcheb(z, 2.26666899049817806459E-8, ref b0, ref b1, ref b2); besselmnextcheb(z, 2.04891858946906374183E-7, ref b0, ref b1, ref b2); besselmnextcheb(z, 2.89137052083475648297E-6, ref b0, ref b1, ref b2); besselmnextcheb(z, 6.88975834691682398426E-5, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.36911647825569408990E-3, ref b0, ref b1, ref b2); besselmnextcheb(z, 8.04490411014108831608E-1, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); result = Math.Exp(x)*v/Math.Sqrt(x); return result; } /************************************************************************* Modified Bessel function of order one Returns modified Bessel function of order one of the argument. The function is defined as i1(x) = -i j1( ix ). The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.9e-15 2.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1985, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besseli1(double x) { double result = 0; double y = 0; double z = 0; double v = 0; double b0 = 0; double b1 = 0; double b2 = 0; z = Math.Abs(x); if( (double)(z)<=(double)(8.0) ) { y = z/2.0-2.0; besselm1firstcheb(2.77791411276104639959E-18, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.11142121435816608115E-17, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.55363195773620046921E-16, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.10559694773538630805E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, 7.60068429473540693410E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, -5.04218550472791168711E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.22379336594557470981E-13, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.98397439776494371520E-12, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.17361862988909016308E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, -6.66348972350202774223E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.62559028155211703701E-10, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.88724975172282928790E-9, ref b0, ref b1, ref b2); besselm1nextcheb(y, 9.38153738649577178388E-9, ref b0, ref b1, ref b2); besselm1nextcheb(y, -4.44505912879632808065E-8, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.00329475355213526229E-7, ref b0, ref b1, ref b2); besselm1nextcheb(y, -8.56872026469545474066E-7, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.47025130813767847674E-6, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.32731636560394358279E-5, ref b0, ref b1, ref b2); besselm1nextcheb(y, 4.78156510755005422638E-5, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.61760815825896745588E-4, ref b0, ref b1, ref b2); besselm1nextcheb(y, 5.12285956168575772895E-4, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.51357245063125314899E-3, ref b0, ref b1, ref b2); besselm1nextcheb(y, 4.15642294431288815669E-3, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.05640848946261981558E-2, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.47264490306265168283E-2, ref b0, ref b1, ref b2); besselm1nextcheb(y, -5.29459812080949914269E-2, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.02643658689847095384E-1, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.76416518357834055153E-1, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.52587186443633654823E-1, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); z = v*z*Math.Exp(z); } else { y = 32.0/z-2.0; besselm1firstcheb(7.51729631084210481353E-18, ref b0, ref b1, ref b2); besselm1nextcheb(y, 4.41434832307170791151E-18, ref b0, ref b1, ref b2); besselm1nextcheb(y, -4.65030536848935832153E-17, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.20952592199342395980E-17, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.96262899764595013876E-16, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.30820231092092828324E-16, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.88035477551078244854E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.81440307243700780478E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.04202769841288027642E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, 4.27244001671195135429E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.10154184277266431302E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, -4.08355111109219731823E-13, ref b0, ref b1, ref b2); besselm1nextcheb(y, -7.19855177624590851209E-13, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.03562854414708950722E-12, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.41258074366137813316E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.25260358301548823856E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.89749581235054123450E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, -5.58974346219658380687E-10, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.83538038596423702205E-9, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.63146884688951950684E-8, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.51223623787020892529E-7, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.88256480887769039346E-6, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.10588938762623716291E-4, ref b0, ref b1, ref b2); besselm1nextcheb(y, -9.76109749136146840777E-3, ref b0, ref b1, ref b2); besselm1nextcheb(y, 7.78576235018280120474E-1, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); z = v*Math.Exp(z)/Math.Sqrt(z); } if( (double)(x)<(double)(0) ) { z = -z; } result = z; return result; } /************************************************************************* Modified Bessel function, second kind, order zero Returns modified Bessel function of the second kind of order zero of the argument. The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Tested at 2000 random points between 0 and 8. Peak absolute error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. Relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.2e-15 1.6e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselk0(double x) { double result = 0; double y = 0; double z = 0; double v = 0; double b0 = 0; double b1 = 0; double b2 = 0; System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselK0: x<=0"); if( (double)(x)<=(double)(2) ) { y = x*x-2.0; besselmfirstcheb(1.37446543561352307156E-16, ref b0, ref b1, ref b2); besselmnextcheb(y, 4.25981614279661018399E-14, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.03496952576338420167E-11, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.90451637722020886025E-9, ref b0, ref b1, ref b2); besselmnextcheb(y, 2.53479107902614945675E-7, ref b0, ref b1, ref b2); besselmnextcheb(y, 2.28621210311945178607E-5, ref b0, ref b1, ref b2); besselmnextcheb(y, 1.26461541144692592338E-3, ref b0, ref b1, ref b2); besselmnextcheb(y, 3.59799365153615016266E-2, ref b0, ref b1, ref b2); besselmnextcheb(y, 3.44289899924628486886E-1, ref b0, ref b1, ref b2); besselmnextcheb(y, -5.35327393233902768720E-1, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); v = v-Math.Log(0.5*x)*besseli0(x); } else { z = 8.0/x-2.0; besselmfirstcheb(5.30043377268626276149E-18, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.64758043015242134646E-17, ref b0, ref b1, ref b2); besselmnextcheb(z, 5.21039150503902756861E-17, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.67823109680541210385E-16, ref b0, ref b1, ref b2); besselmnextcheb(z, 5.51205597852431940784E-16, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.84859337734377901440E-15, ref b0, ref b1, ref b2); besselmnextcheb(z, 6.34007647740507060557E-15, ref b0, ref b1, ref b2); besselmnextcheb(z, -2.22751332699166985548E-14, ref b0, ref b1, ref b2); besselmnextcheb(z, 8.03289077536357521100E-14, ref b0, ref b1, ref b2); besselmnextcheb(z, -2.98009692317273043925E-13, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.14034058820847496303E-12, ref b0, ref b1, ref b2); besselmnextcheb(z, -4.51459788337394416547E-12, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.85594911495471785253E-11, ref b0, ref b1, ref b2); besselmnextcheb(z, -7.95748924447710747776E-11, ref b0, ref b1, ref b2); besselmnextcheb(z, 3.57739728140030116597E-10, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.69753450938905987466E-9, ref b0, ref b1, ref b2); besselmnextcheb(z, 8.57403401741422608519E-9, ref b0, ref b1, ref b2); besselmnextcheb(z, -4.66048989768794782956E-8, ref b0, ref b1, ref b2); besselmnextcheb(z, 2.76681363944501510342E-7, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.83175552271911948767E-6, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.39498137188764993662E-5, ref b0, ref b1, ref b2); besselmnextcheb(z, -1.28495495816278026384E-4, ref b0, ref b1, ref b2); besselmnextcheb(z, 1.56988388573005337491E-3, ref b0, ref b1, ref b2); besselmnextcheb(z, -3.14481013119645005427E-2, ref b0, ref b1, ref b2); besselmnextcheb(z, 2.44030308206595545468E0, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); v = v*Math.Exp(-x)/Math.Sqrt(x); } result = v; return result; } /************************************************************************* Modified Bessel function, second kind, order one Computes the modified Bessel function of the second kind of order one of the argument. The range is partitioned into the two intervals [0,2] and (2, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.2e-15 1.6e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselk1(double x) { double result = 0; double y = 0; double z = 0; double v = 0; double b0 = 0; double b1 = 0; double b2 = 0; z = 0.5*x; System.Diagnostics.Debug.Assert((double)(z)>(double)(0), "Domain error in K1"); if( (double)(x)<=(double)(2) ) { y = x*x-2.0; besselm1firstcheb(-7.02386347938628759343E-18, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.42744985051936593393E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, -6.66690169419932900609E-13, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.41148839263352776110E-10, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.21338763073472585583E-8, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.43340614156596823496E-6, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.73028895751305206302E-4, ref b0, ref b1, ref b2); besselm1nextcheb(y, -6.97572385963986435018E-3, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.22611180822657148235E-1, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.53155960776544875667E-1, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.52530022733894777053E0, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); result = Math.Log(z)*besseli1(x)+v/x; } else { y = 8.0/x-2.0; besselm1firstcheb(-5.75674448366501715755E-18, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.79405087314755922667E-17, ref b0, ref b1, ref b2); besselm1nextcheb(y, -5.68946255844285935196E-17, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.83809354436663880070E-16, ref b0, ref b1, ref b2); besselm1nextcheb(y, -6.05704724837331885336E-16, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.03870316562433424052E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, -7.01983709041831346144E-15, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.47715442448130437068E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, -8.97670518232499435011E-14, ref b0, ref b1, ref b2); besselm1nextcheb(y, 3.34841966607842919884E-13, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.28917396095102890680E-12, ref b0, ref b1, ref b2); besselm1nextcheb(y, 5.13963967348173025100E-12, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.12996783842756842877E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, 9.21831518760500529508E-11, ref b0, ref b1, ref b2); besselm1nextcheb(y, -4.19035475934189648750E-10, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.01504975519703286596E-9, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.03457624656780970260E-8, ref b0, ref b1, ref b2); besselm1nextcheb(y, 5.74108412545004946722E-8, ref b0, ref b1, ref b2); besselm1nextcheb(y, -3.50196060308781257119E-7, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.40648494783721712015E-6, ref b0, ref b1, ref b2); besselm1nextcheb(y, -1.93619797416608296024E-5, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.95215518471351631108E-4, ref b0, ref b1, ref b2); besselm1nextcheb(y, -2.85781685962277938680E-3, ref b0, ref b1, ref b2); besselm1nextcheb(y, 1.03923736576817238437E-1, ref b0, ref b1, ref b2); besselm1nextcheb(y, 2.72062619048444266945E0, ref b0, ref b1, ref b2); v = 0.5*(b0-b2); result = Math.Exp(-x)*v/Math.Sqrt(x); } return result; } /************************************************************************* Modified Bessel function, second kind, integer order Returns modified Bessel function of the second kind of order n of the argument. The range is partitioned into the two intervals [0,9.55] and (9.55, infinity). An ascending power series is used in the low range, and an asymptotic expansion in the high range. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,30 90000 1.8e-8 3.0e-10 Error is high only near the crossover point x = 9.55 between the two expansions used. Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier *************************************************************************/ public static double besselkn(int nn, double x) { double result = 0; double k = 0; double kf = 0; double nk1f = 0; double nkf = 0; double zn = 0; double t = 0; double s = 0; double z0 = 0; double z = 0; double ans = 0; double fn = 0; double pn = 0; double pk = 0; double zmn = 0; double tlg = 0; double tox = 0; int i = 0; int n = 0; double eul = 0; eul = 5.772156649015328606065e-1; if( nn<0 ) { n = -nn; } else { n = nn; } System.Diagnostics.Debug.Assert(n<=31, "Overflow in BesselKN"); System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselKN"); if( (double)(x)<=(double)(9.55) ) { ans = 0.0; z0 = 0.25*x*x; fn = 1.0; pn = 0.0; zmn = 1.0; tox = 2.0/x; if( n>0 ) { pn = -eul; k = 1.0; for(i=1; i<=n-1; i++) { pn = pn+1.0/k; k = k+1.0; fn = fn*k; } zmn = tox; if( n==1 ) { ans = 1.0/x; } else { nk1f = fn/n; kf = 1.0; s = nk1f; z = -z0; zn = 1.0; for(i=1; i<=n-1; i++) { nk1f = nk1f/(n-i); kf = kf*i; zn = zn*z; t = nk1f*zn/kf; s = s+t; System.Diagnostics.Debug.Assert((double)(AP.Math.MaxRealNumber-Math.Abs(t))>(double)(Math.Abs(s)), "Overflow in BesselKN"); System.Diagnostics.Debug.Assert(!((double)(tox)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/tox)<(double)(zmn)), "Overflow in BesselKN"); zmn = zmn*tox; } s = s*0.5; t = Math.Abs(s); System.Diagnostics.Debug.Assert(!((double)(zmn)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/zmn)<(double)(t)), "Overflow in BesselKN"); System.Diagnostics.Debug.Assert(!((double)(t)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/t)<(double)(zmn)), "Overflow in BesselKN"); ans = s*zmn; } } tlg = 2.0*Math.Log(0.5*x); pk = -eul; if( n==0 ) { pn = pk; t = 1.0; } else { pn = pn+1.0/n; t = 1.0/fn; } s = (pk+pn-tlg)*t; k = 1.0; do { t = t*(z0/(k*(k+n))); pk = pk+1.0/k; pn = pn+1.0/(k+n); s = s+(pk+pn-tlg)*t; k = k+1.0; } while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) ); s = 0.5*s/zmn; if( n%2!=0 ) { s = -s; } ans = ans+s; result = ans; return result; } if( (double)(x)>(double)(Math.Log(AP.Math.MaxRealNumber)) ) { result = 0; return result; } k = n; pn = 4.0*k*k; pk = 1.0; z0 = 8.0*x; fn = 1.0; t = 1.0; s = t; nkf = AP.Math.MaxRealNumber; i = 0; do { z = pn-pk*pk; t = t*z/(fn*z0); nk1f = Math.Abs(t); if( i>=n & (double)(nk1f)>(double)(nkf) ) { break; } nkf = nk1f; s = s+t; fn = fn+1.0; pk = pk+2.0; i = i+1; } while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) ); result = Math.Exp(-x)*Math.Sqrt(Math.PI/(2.0*x))*s; return result; } /************************************************************************* Internal subroutine Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ private static void besselmfirstcheb(double c, ref double b0, ref double b1, ref double b2) { b0 = c; b1 = 0.0; b2 = 0.0; } /************************************************************************* Internal subroutine Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ private static void besselmnextcheb(double x, double c, ref double b0, ref double b1, ref double b2) { b2 = b1; b1 = b0; b0 = x*b1-b2+c; } /************************************************************************* Internal subroutine Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ private static void besselm1firstcheb(double c, ref double b0, ref double b1, ref double b2) { b0 = c; b1 = 0.0; b2 = 0.0; } /************************************************************************* Internal subroutine Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ private static void besselm1nextcheb(double x, double c, ref double b0, ref double b1, ref double b2) { b2 = b1; b1 = b0; b0 = x*b1-b2+c; } private static void besselasympt0(double x, ref double pzero, ref double qzero) { double xsq = 0; double p2 = 0; double q2 = 0; double p3 = 0; double q3 = 0; xsq = 64.0/(x*x); p2 = 0.0; p2 = 2485.271928957404011288128951+xsq*p2; p2 = 153982.6532623911470917825993+xsq*p2; p2 = 2016135.283049983642487182349+xsq*p2; p2 = 8413041.456550439208464315611+xsq*p2; p2 = 12332384.76817638145232406055+xsq*p2; p2 = 5393485.083869438325262122897+xsq*p2; q2 = 1.0; q2 = 2615.700736920839685159081813+xsq*q2; q2 = 156001.7276940030940592769933+xsq*q2; q2 = 2025066.801570134013891035236+xsq*q2; q2 = 8426449.050629797331554404810+xsq*q2; q2 = 12338310.22786324960844856182+xsq*q2; q2 = 5393485.083869438325560444960+xsq*q2; p3 = -0.0; p3 = -4.887199395841261531199129300+xsq*p3; p3 = -226.2630641933704113967255053+xsq*p3; p3 = -2365.956170779108192723612816+xsq*p3; p3 = -8239.066313485606568803548860+xsq*p3; p3 = -10381.41698748464093880530341+xsq*p3; p3 = -3984.617357595222463506790588+xsq*p3; q3 = 1.0; q3 = 408.7714673983499223402830260+xsq*q3; q3 = 15704.89191515395519392882766+xsq*q3; q3 = 156021.3206679291652539287109+xsq*q3; q3 = 533291.3634216897168722255057+xsq*q3; q3 = 666745.4239319826986004038103+xsq*q3; q3 = 255015.5108860942382983170882+xsq*q3; pzero = p2/q2; qzero = 8*p3/q3/x; } private static void besselasympt1(double x, ref double pzero, ref double qzero) { double xsq = 0; double p2 = 0; double q2 = 0; double p3 = 0; double q3 = 0; xsq = 64.0/(x*x); p2 = -1611.616644324610116477412898; p2 = -109824.0554345934672737413139+xsq*p2; p2 = -1523529.351181137383255105722+xsq*p2; p2 = -6603373.248364939109255245434+xsq*p2; p2 = -9942246.505077641195658377899+xsq*p2; p2 = -4435757.816794127857114720794+xsq*p2; q2 = 1.0; q2 = -1455.009440190496182453565068+xsq*q2; q2 = -107263.8599110382011903063867+xsq*q2; q2 = -1511809.506634160881644546358+xsq*q2; q2 = -6585339.479723087072826915069+xsq*q2; q2 = -9934124.389934585658967556309+xsq*q2; q2 = -4435757.816794127856828016962+xsq*q2; p3 = 35.26513384663603218592175580; p3 = 1706.375429020768002061283546+xsq*p3; p3 = 18494.26287322386679652009819+xsq*p3; p3 = 66178.83658127083517939992166+xsq*p3; p3 = 85145.16067533570196555001171+xsq*p3; p3 = 33220.91340985722351859704442+xsq*p3; q3 = 1.0; q3 = 863.8367769604990967475517183+xsq*q3; q3 = 37890.22974577220264142952256+xsq*q3; q3 = 400294.4358226697511708610813+xsq*q3; q3 = 1419460.669603720892855755253+xsq*q3; q3 = 1819458.042243997298924553839+xsq*q3; q3 = 708712.8194102874357377502472+xsq*q3; pzero = p2/q2; qzero = 8*p3/q3/x; } } }