Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/bessel.cs @ 17712

Last change on this file since 17712 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 50.8 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace alglib
31{
32    public class bessel
33    {
34        /*************************************************************************
35        Bessel function of order zero
36
37        Returns Bessel function of order zero of the argument.
38
39        The domain is divided into the intervals [0, 5] and
40        (5, infinity). In the first interval the following rational
41        approximation is used:
42
43
44               2         2
45        (w - r  ) (w - r  ) P (w) / Q (w)
46              1         2    3       8
47
48                   2
49        where w = x  and the two r's are zeros of the function.
50
51        In the second interval, the Hankel asymptotic expansion
52        is employed with two rational functions of degree 6/6
53        and 7/7.
54
55        ACCURACY:
56
57                             Absolute error:
58        arithmetic   domain     # trials      peak         rms
59           IEEE      0, 30       60000       4.2e-16     1.1e-16
60
61        Cephes Math Library Release 2.8:  June, 2000
62        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
63        *************************************************************************/
64        public static double besselj0(double x)
65        {
66            double result = 0;
67            double xsq = 0;
68            double nn = 0;
69            double pzero = 0;
70            double qzero = 0;
71            double p1 = 0;
72            double q1 = 0;
73
74            if( (double)(x)<(double)(0) )
75            {
76                x = -x;
77            }
78            if( (double)(x)>(double)(8.0) )
79            {
80                besselasympt0(x, ref pzero, ref qzero);
81                nn = x-Math.PI/4;
82                result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn));
83                return result;
84            }
85            xsq = AP.Math.Sqr(x);
86            p1 = 26857.86856980014981415848441;
87            p1 = -40504123.71833132706360663322+xsq*p1;
88            p1 = 25071582855.36881945555156435+xsq*p1;
89            p1 = -8085222034853.793871199468171+xsq*p1;
90            p1 = 1434354939140344.111664316553+xsq*p1;
91            p1 = -136762035308817138.6865416609+xsq*p1;
92            p1 = 6382059341072356562.289432465+xsq*p1;
93            p1 = -117915762910761053603.8440800+xsq*p1;
94            p1 = 493378725179413356181.6813446+xsq*p1;
95            q1 = 1.0;
96            q1 = 1363.063652328970604442810507+xsq*q1;
97            q1 = 1114636.098462985378182402543+xsq*q1;
98            q1 = 669998767.2982239671814028660+xsq*q1;
99            q1 = 312304311494.1213172572469442+xsq*q1;
100            q1 = 112775673967979.8507056031594+xsq*q1;
101            q1 = 30246356167094626.98627330784+xsq*q1;
102            q1 = 5428918384092285160.200195092+xsq*q1;
103            q1 = 493378725179413356211.3278438+xsq*q1;
104            result = p1/q1;
105            return result;
106        }
107
108
109        /*************************************************************************
110        Bessel function of order one
111
112        Returns Bessel function of order one of the argument.
113
114        The domain is divided into the intervals [0, 8] and
115        (8, infinity). In the first interval a 24 term Chebyshev
116        expansion is used. In the second, the asymptotic
117        trigonometric representation is employed using two
118        rational functions of degree 5/5.
119
120        ACCURACY:
121
122                             Absolute error:
123        arithmetic   domain      # trials      peak         rms
124           IEEE      0, 30       30000       2.6e-16     1.1e-16
125
126        Cephes Math Library Release 2.8:  June, 2000
127        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
128        *************************************************************************/
129        public static double besselj1(double x)
130        {
131            double result = 0;
132            double s = 0;
133            double xsq = 0;
134            double nn = 0;
135            double pzero = 0;
136            double qzero = 0;
137            double p1 = 0;
138            double q1 = 0;
139
140            s = Math.Sign(x);
141            if( (double)(x)<(double)(0) )
142            {
143                x = -x;
144            }
145            if( (double)(x)>(double)(8.0) )
146            {
147                besselasympt1(x, ref pzero, ref qzero);
148                nn = x-3*Math.PI/4;
149                result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn));
150                if( (double)(s)<(double)(0) )
151                {
152                    result = -result;
153                }
154                return result;
155            }
156            xsq = AP.Math.Sqr(x);
157            p1 = 2701.122710892323414856790990;
158            p1 = -4695753.530642995859767162166+xsq*p1;
159            p1 = 3413234182.301700539091292655+xsq*p1;
160            p1 = -1322983480332.126453125473247+xsq*p1;
161            p1 = 290879526383477.5409737601689+xsq*p1;
162            p1 = -35888175699101060.50743641413+xsq*p1;
163            p1 = 2316433580634002297.931815435+xsq*p1;
164            p1 = -66721065689249162980.20941484+xsq*p1;
165            p1 = 581199354001606143928.050809+xsq*p1;
166            q1 = 1.0;
167            q1 = 1606.931573481487801970916749+xsq*q1;
168            q1 = 1501793.594998585505921097578+xsq*q1;
169            q1 = 1013863514.358673989967045588+xsq*q1;
170            q1 = 524371026216.7649715406728642+xsq*q1;
171            q1 = 208166122130760.7351240184229+xsq*q1;
172            q1 = 60920613989175217.46105196863+xsq*q1;
173            q1 = 11857707121903209998.37113348+xsq*q1;
174            q1 = 1162398708003212287858.529400+xsq*q1;
175            result = s*x*p1/q1;
176            return result;
177        }
178
179
180        /*************************************************************************
181        Bessel function of integer order
182
183        Returns Bessel function of order n, where n is a
184        (possibly negative) integer.
185
186        The ratio of jn(x) to j0(x) is computed by backward
187        recurrence.  First the ratio jn/jn-1 is found by a
188        continued fraction expansion.  Then the recurrence
189        relating successive orders is applied until j0 or j1 is
190        reached.
191
192        If n = 0 or 1 the routine for j0 or j1 is called
193        directly.
194
195        ACCURACY:
196
197                             Absolute error:
198        arithmetic   range      # trials      peak         rms
199           IEEE      0, 30        5000       4.4e-16     7.9e-17
200
201
202        Not suitable for large n or x. Use jv() (fractional order) instead.
203
204        Cephes Math Library Release 2.8:  June, 2000
205        Copyright 1984, 1987, 2000 by Stephen L. Moshier
206        *************************************************************************/
207        public static double besseljn(int n,
208            double x)
209        {
210            double result = 0;
211            double pkm2 = 0;
212            double pkm1 = 0;
213            double pk = 0;
214            double xk = 0;
215            double r = 0;
216            double ans = 0;
217            int k = 0;
218            int sg = 0;
219
220            if( n<0 )
221            {
222                n = -n;
223                if( n%2==0 )
224                {
225                    sg = 1;
226                }
227                else
228                {
229                    sg = -1;
230                }
231            }
232            else
233            {
234                sg = 1;
235            }
236            if( (double)(x)<(double)(0) )
237            {
238                if( n%2!=0 )
239                {
240                    sg = -sg;
241                }
242                x = -x;
243            }
244            if( n==0 )
245            {
246                result = sg*besselj0(x);
247                return result;
248            }
249            if( n==1 )
250            {
251                result = sg*besselj1(x);
252                return result;
253            }
254            if( n==2 )
255            {
256                if( (double)(x)==(double)(0) )
257                {
258                    result = 0;
259                }
260                else
261                {
262                    result = sg*(2.0*besselj1(x)/x-besselj0(x));
263                }
264                return result;
265            }
266            if( (double)(x)<(double)(AP.Math.MachineEpsilon) )
267            {
268                result = 0;
269                return result;
270            }
271            k = 53;
272            pk = 2*(n+k);
273            ans = pk;
274            xk = x*x;
275            do
276            {
277                pk = pk-2.0;
278                ans = pk-xk/ans;
279                k = k-1;
280            }
281            while( k!=0 );
282            ans = x/ans;
283            pk = 1.0;
284            pkm1 = 1.0/ans;
285            k = n-1;
286            r = 2*k;
287            do
288            {
289                pkm2 = (pkm1*r-pk*x)/x;
290                pk = pkm1;
291                pkm1 = pkm2;
292                r = r-2.0;
293                k = k-1;
294            }
295            while( k!=0 );
296            if( (double)(Math.Abs(pk))>(double)(Math.Abs(pkm1)) )
297            {
298                ans = besselj1(x)/pk;
299            }
300            else
301            {
302                ans = besselj0(x)/pkm1;
303            }
304            result = sg*ans;
305            return result;
306        }
307
308
309        /*************************************************************************
310        Bessel function of the second kind, order zero
311
312        Returns Bessel function of the second kind, of order
313        zero, of the argument.
314
315        The domain is divided into the intervals [0, 5] and
316        (5, infinity). In the first interval a rational approximation
317        R(x) is employed to compute
318          y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
319        Thus a call to j0() is required.
320
321        In the second interval, the Hankel asymptotic expansion
322        is employed with two rational functions of degree 6/6
323        and 7/7.
324
325
326
327        ACCURACY:
328
329         Absolute error, when y0(x) < 1; else relative error:
330
331        arithmetic   domain     # trials      peak         rms
332           IEEE      0, 30       30000       1.3e-15     1.6e-16
333
334        Cephes Math Library Release 2.8:  June, 2000
335        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
336        *************************************************************************/
337        public static double bessely0(double x)
338        {
339            double result = 0;
340            double nn = 0;
341            double xsq = 0;
342            double pzero = 0;
343            double qzero = 0;
344            double p4 = 0;
345            double q4 = 0;
346
347            if( (double)(x)>(double)(8.0) )
348            {
349                besselasympt0(x, ref pzero, ref qzero);
350                nn = x-Math.PI/4;
351                result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn));
352                return result;
353            }
354            xsq = AP.Math.Sqr(x);
355            p4 = -41370.35497933148554125235152;
356            p4 = 59152134.65686889654273830069+xsq*p4;
357            p4 = -34363712229.79040378171030138+xsq*p4;
358            p4 = 10255208596863.94284509167421+xsq*p4;
359            p4 = -1648605817185729.473122082537+xsq*p4;
360            p4 = 137562431639934407.8571335453+xsq*p4;
361            p4 = -5247065581112764941.297350814+xsq*p4;
362            p4 = 65874732757195549259.99402049+xsq*p4;
363            p4 = -27502866786291095837.01933175+xsq*p4;
364            q4 = 1.0;
365            q4 = 1282.452772478993804176329391+xsq*q4;
366            q4 = 1001702.641288906265666651753+xsq*q4;
367            q4 = 579512264.0700729537480087915+xsq*q4;
368            q4 = 261306575504.1081249568482092+xsq*q4;
369            q4 = 91620380340751.85262489147968+xsq*q4;
370            q4 = 23928830434997818.57439356652+xsq*q4;
371            q4 = 4192417043410839973.904769661+xsq*q4;
372            q4 = 372645883898616588198.9980+xsq*q4;
373            result = p4/q4+2/Math.PI*besselj0(x)*Math.Log(x);
374            return result;
375        }
376
377
378        /*************************************************************************
379        Bessel function of second kind of order one
380
381        Returns Bessel function of the second kind of order one
382        of the argument.
383
384        The domain is divided into the intervals [0, 8] and
385        (8, infinity). In the first interval a 25 term Chebyshev
386        expansion is used, and a call to j1() is required.
387        In the second, the asymptotic trigonometric representation
388        is employed using two rational functions of degree 5/5.
389
390        ACCURACY:
391
392                             Absolute error:
393        arithmetic   domain      # trials      peak         rms
394           IEEE      0, 30       30000       1.0e-15     1.3e-16
395
396        Cephes Math Library Release 2.8:  June, 2000
397        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
398        *************************************************************************/
399        public static double bessely1(double x)
400        {
401            double result = 0;
402            double nn = 0;
403            double xsq = 0;
404            double pzero = 0;
405            double qzero = 0;
406            double p4 = 0;
407            double q4 = 0;
408
409            if( (double)(x)>(double)(8.0) )
410            {
411                besselasympt1(x, ref pzero, ref qzero);
412                nn = x-3*Math.PI/4;
413                result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn));
414                return result;
415            }
416            xsq = AP.Math.Sqr(x);
417            p4 = -2108847.540133123652824139923;
418            p4 = 3639488548.124002058278999428+xsq*p4;
419            p4 = -2580681702194.450950541426399+xsq*p4;
420            p4 = 956993023992168.3481121552788+xsq*p4;
421            p4 = -196588746272214065.8820322248+xsq*p4;
422            p4 = 21931073399177975921.11427556+xsq*p4;
423            p4 = -1212297555414509577913.561535+xsq*p4;
424            p4 = 26554738314348543268942.48968+xsq*p4;
425            p4 = -99637534243069222259967.44354+xsq*p4;
426            q4 = 1.0;
427            q4 = 1612.361029677000859332072312+xsq*q4;
428            q4 = 1563282.754899580604737366452+xsq*q4;
429            q4 = 1128686837.169442121732366891+xsq*q4;
430            q4 = 646534088126.5275571961681500+xsq*q4;
431            q4 = 297663212564727.6729292742282+xsq*q4;
432            q4 = 108225825940881955.2553850180+xsq*q4;
433            q4 = 29549879358971486742.90758119+xsq*q4;
434            q4 = 5435310377188854170800.653097+xsq*q4;
435            q4 = 508206736694124324531442.4152+xsq*q4;
436            result = x*p4/q4+2/Math.PI*(besselj1(x)*Math.Log(x)-1/x);
437            return result;
438        }
439
440
441        /*************************************************************************
442        Bessel function of second kind of integer order
443
444        Returns Bessel function of order n, where n is a
445        (possibly negative) integer.
446
447        The function is evaluated by forward recurrence on
448        n, starting with values computed by the routines
449        y0() and y1().
450
451        If n = 0 or 1 the routine for y0 or y1 is called
452        directly.
453
454        ACCURACY:
455                             Absolute error, except relative
456                             when y > 1:
457        arithmetic   domain     # trials      peak         rms
458           IEEE      0, 30       30000       3.4e-15     4.3e-16
459
460        Cephes Math Library Release 2.8:  June, 2000
461        Copyright 1984, 1987, 2000 by Stephen L. Moshier
462        *************************************************************************/
463        public static double besselyn(int n,
464            double x)
465        {
466            double result = 0;
467            int i = 0;
468            double a = 0;
469            double b = 0;
470            double tmp = 0;
471            double s = 0;
472
473            s = 1;
474            if( n<0 )
475            {
476                n = -n;
477                if( n%2!=0 )
478                {
479                    s = -1;
480                }
481            }
482            if( n==0 )
483            {
484                result = bessely0(x);
485                return result;
486            }
487            if( n==1 )
488            {
489                result = s*bessely1(x);
490                return result;
491            }
492            a = bessely0(x);
493            b = bessely1(x);
494            for(i=1; i<=n-1; i++)
495            {
496                tmp = b;
497                b = 2*i/x*b-a;
498                a = tmp;
499            }
500            result = s*b;
501            return result;
502        }
503
504
505        /*************************************************************************
506        Modified Bessel function of order zero
507
508        Returns modified Bessel function of order zero of the
509        argument.
510
511        The function is defined as i0(x) = j0( ix ).
512
513        The range is partitioned into the two intervals [0,8] and
514        (8, infinity).  Chebyshev polynomial expansions are employed
515        in each interval.
516
517        ACCURACY:
518
519                             Relative error:
520        arithmetic   domain     # trials      peak         rms
521           IEEE      0,30        30000       5.8e-16     1.4e-16
522
523        Cephes Math Library Release 2.8:  June, 2000
524        Copyright 1984, 1987, 2000 by Stephen L. Moshier
525        *************************************************************************/
526        public static double besseli0(double x)
527        {
528            double result = 0;
529            double y = 0;
530            double v = 0;
531            double z = 0;
532            double b0 = 0;
533            double b1 = 0;
534            double b2 = 0;
535
536            if( (double)(x)<(double)(0) )
537            {
538                x = -x;
539            }
540            if( (double)(x)<=(double)(8.0) )
541            {
542                y = x/2.0-2.0;
543                besselmfirstcheb(-4.41534164647933937950E-18, ref b0, ref b1, ref b2);
544                besselmnextcheb(y, 3.33079451882223809783E-17, ref b0, ref b1, ref b2);
545                besselmnextcheb(y, -2.43127984654795469359E-16, ref b0, ref b1, ref b2);
546                besselmnextcheb(y, 1.71539128555513303061E-15, ref b0, ref b1, ref b2);
547                besselmnextcheb(y, -1.16853328779934516808E-14, ref b0, ref b1, ref b2);
548                besselmnextcheb(y, 7.67618549860493561688E-14, ref b0, ref b1, ref b2);
549                besselmnextcheb(y, -4.85644678311192946090E-13, ref b0, ref b1, ref b2);
550                besselmnextcheb(y, 2.95505266312963983461E-12, ref b0, ref b1, ref b2);
551                besselmnextcheb(y, -1.72682629144155570723E-11, ref b0, ref b1, ref b2);
552                besselmnextcheb(y, 9.67580903537323691224E-11, ref b0, ref b1, ref b2);
553                besselmnextcheb(y, -5.18979560163526290666E-10, ref b0, ref b1, ref b2);
554                besselmnextcheb(y, 2.65982372468238665035E-9, ref b0, ref b1, ref b2);
555                besselmnextcheb(y, -1.30002500998624804212E-8, ref b0, ref b1, ref b2);
556                besselmnextcheb(y, 6.04699502254191894932E-8, ref b0, ref b1, ref b2);
557                besselmnextcheb(y, -2.67079385394061173391E-7, ref b0, ref b1, ref b2);
558                besselmnextcheb(y, 1.11738753912010371815E-6, ref b0, ref b1, ref b2);
559                besselmnextcheb(y, -4.41673835845875056359E-6, ref b0, ref b1, ref b2);
560                besselmnextcheb(y, 1.64484480707288970893E-5, ref b0, ref b1, ref b2);
561                besselmnextcheb(y, -5.75419501008210370398E-5, ref b0, ref b1, ref b2);
562                besselmnextcheb(y, 1.88502885095841655729E-4, ref b0, ref b1, ref b2);
563                besselmnextcheb(y, -5.76375574538582365885E-4, ref b0, ref b1, ref b2);
564                besselmnextcheb(y, 1.63947561694133579842E-3, ref b0, ref b1, ref b2);
565                besselmnextcheb(y, -4.32430999505057594430E-3, ref b0, ref b1, ref b2);
566                besselmnextcheb(y, 1.05464603945949983183E-2, ref b0, ref b1, ref b2);
567                besselmnextcheb(y, -2.37374148058994688156E-2, ref b0, ref b1, ref b2);
568                besselmnextcheb(y, 4.93052842396707084878E-2, ref b0, ref b1, ref b2);
569                besselmnextcheb(y, -9.49010970480476444210E-2, ref b0, ref b1, ref b2);
570                besselmnextcheb(y, 1.71620901522208775349E-1, ref b0, ref b1, ref b2);
571                besselmnextcheb(y, -3.04682672343198398683E-1, ref b0, ref b1, ref b2);
572                besselmnextcheb(y, 6.76795274409476084995E-1, ref b0, ref b1, ref b2);
573                v = 0.5*(b0-b2);
574                result = Math.Exp(x)*v;
575                return result;
576            }
577            z = 32.0/x-2.0;
578            besselmfirstcheb(-7.23318048787475395456E-18, ref b0, ref b1, ref b2);
579            besselmnextcheb(z, -4.83050448594418207126E-18, ref b0, ref b1, ref b2);
580            besselmnextcheb(z, 4.46562142029675999901E-17, ref b0, ref b1, ref b2);
581            besselmnextcheb(z, 3.46122286769746109310E-17, ref b0, ref b1, ref b2);
582            besselmnextcheb(z, -2.82762398051658348494E-16, ref b0, ref b1, ref b2);
583            besselmnextcheb(z, -3.42548561967721913462E-16, ref b0, ref b1, ref b2);
584            besselmnextcheb(z, 1.77256013305652638360E-15, ref b0, ref b1, ref b2);
585            besselmnextcheb(z, 3.81168066935262242075E-15, ref b0, ref b1, ref b2);
586            besselmnextcheb(z, -9.55484669882830764870E-15, ref b0, ref b1, ref b2);
587            besselmnextcheb(z, -4.15056934728722208663E-14, ref b0, ref b1, ref b2);
588            besselmnextcheb(z, 1.54008621752140982691E-14, ref b0, ref b1, ref b2);
589            besselmnextcheb(z, 3.85277838274214270114E-13, ref b0, ref b1, ref b2);
590            besselmnextcheb(z, 7.18012445138366623367E-13, ref b0, ref b1, ref b2);
591            besselmnextcheb(z, -1.79417853150680611778E-12, ref b0, ref b1, ref b2);
592            besselmnextcheb(z, -1.32158118404477131188E-11, ref b0, ref b1, ref b2);
593            besselmnextcheb(z, -3.14991652796324136454E-11, ref b0, ref b1, ref b2);
594            besselmnextcheb(z, 1.18891471078464383424E-11, ref b0, ref b1, ref b2);
595            besselmnextcheb(z, 4.94060238822496958910E-10, ref b0, ref b1, ref b2);
596            besselmnextcheb(z, 3.39623202570838634515E-9, ref b0, ref b1, ref b2);
597            besselmnextcheb(z, 2.26666899049817806459E-8, ref b0, ref b1, ref b2);
598            besselmnextcheb(z, 2.04891858946906374183E-7, ref b0, ref b1, ref b2);
599            besselmnextcheb(z, 2.89137052083475648297E-6, ref b0, ref b1, ref b2);
600            besselmnextcheb(z, 6.88975834691682398426E-5, ref b0, ref b1, ref b2);
601            besselmnextcheb(z, 3.36911647825569408990E-3, ref b0, ref b1, ref b2);
602            besselmnextcheb(z, 8.04490411014108831608E-1, ref b0, ref b1, ref b2);
603            v = 0.5*(b0-b2);
604            result = Math.Exp(x)*v/Math.Sqrt(x);
605            return result;
606        }
607
608
609        /*************************************************************************
610        Modified Bessel function of order one
611
612        Returns modified Bessel function of order one of the
613        argument.
614
615        The function is defined as i1(x) = -i j1( ix ).
616
617        The range is partitioned into the two intervals [0,8] and
618        (8, infinity).  Chebyshev polynomial expansions are employed
619        in each interval.
620
621        ACCURACY:
622
623                             Relative error:
624        arithmetic   domain     # trials      peak         rms
625           IEEE      0, 30       30000       1.9e-15     2.1e-16
626
627        Cephes Math Library Release 2.8:  June, 2000
628        Copyright 1985, 1987, 2000 by Stephen L. Moshier
629        *************************************************************************/
630        public static double besseli1(double x)
631        {
632            double result = 0;
633            double y = 0;
634            double z = 0;
635            double v = 0;
636            double b0 = 0;
637            double b1 = 0;
638            double b2 = 0;
639
640            z = Math.Abs(x);
641            if( (double)(z)<=(double)(8.0) )
642            {
643                y = z/2.0-2.0;
644                besselm1firstcheb(2.77791411276104639959E-18, ref b0, ref b1, ref b2);
645                besselm1nextcheb(y, -2.11142121435816608115E-17, ref b0, ref b1, ref b2);
646                besselm1nextcheb(y, 1.55363195773620046921E-16, ref b0, ref b1, ref b2);
647                besselm1nextcheb(y, -1.10559694773538630805E-15, ref b0, ref b1, ref b2);
648                besselm1nextcheb(y, 7.60068429473540693410E-15, ref b0, ref b1, ref b2);
649                besselm1nextcheb(y, -5.04218550472791168711E-14, ref b0, ref b1, ref b2);
650                besselm1nextcheb(y, 3.22379336594557470981E-13, ref b0, ref b1, ref b2);
651                besselm1nextcheb(y, -1.98397439776494371520E-12, ref b0, ref b1, ref b2);
652                besselm1nextcheb(y, 1.17361862988909016308E-11, ref b0, ref b1, ref b2);
653                besselm1nextcheb(y, -6.66348972350202774223E-11, ref b0, ref b1, ref b2);
654                besselm1nextcheb(y, 3.62559028155211703701E-10, ref b0, ref b1, ref b2);
655                besselm1nextcheb(y, -1.88724975172282928790E-9, ref b0, ref b1, ref b2);
656                besselm1nextcheb(y, 9.38153738649577178388E-9, ref b0, ref b1, ref b2);
657                besselm1nextcheb(y, -4.44505912879632808065E-8, ref b0, ref b1, ref b2);
658                besselm1nextcheb(y, 2.00329475355213526229E-7, ref b0, ref b1, ref b2);
659                besselm1nextcheb(y, -8.56872026469545474066E-7, ref b0, ref b1, ref b2);
660                besselm1nextcheb(y, 3.47025130813767847674E-6, ref b0, ref b1, ref b2);
661                besselm1nextcheb(y, -1.32731636560394358279E-5, ref b0, ref b1, ref b2);
662                besselm1nextcheb(y, 4.78156510755005422638E-5, ref b0, ref b1, ref b2);
663                besselm1nextcheb(y, -1.61760815825896745588E-4, ref b0, ref b1, ref b2);
664                besselm1nextcheb(y, 5.12285956168575772895E-4, ref b0, ref b1, ref b2);
665                besselm1nextcheb(y, -1.51357245063125314899E-3, ref b0, ref b1, ref b2);
666                besselm1nextcheb(y, 4.15642294431288815669E-3, ref b0, ref b1, ref b2);
667                besselm1nextcheb(y, -1.05640848946261981558E-2, ref b0, ref b1, ref b2);
668                besselm1nextcheb(y, 2.47264490306265168283E-2, ref b0, ref b1, ref b2);
669                besselm1nextcheb(y, -5.29459812080949914269E-2, ref b0, ref b1, ref b2);
670                besselm1nextcheb(y, 1.02643658689847095384E-1, ref b0, ref b1, ref b2);
671                besselm1nextcheb(y, -1.76416518357834055153E-1, ref b0, ref b1, ref b2);
672                besselm1nextcheb(y, 2.52587186443633654823E-1, ref b0, ref b1, ref b2);
673                v = 0.5*(b0-b2);
674                z = v*z*Math.Exp(z);
675            }
676            else
677            {
678                y = 32.0/z-2.0;
679                besselm1firstcheb(7.51729631084210481353E-18, ref b0, ref b1, ref b2);
680                besselm1nextcheb(y, 4.41434832307170791151E-18, ref b0, ref b1, ref b2);
681                besselm1nextcheb(y, -4.65030536848935832153E-17, ref b0, ref b1, ref b2);
682                besselm1nextcheb(y, -3.20952592199342395980E-17, ref b0, ref b1, ref b2);
683                besselm1nextcheb(y, 2.96262899764595013876E-16, ref b0, ref b1, ref b2);
684                besselm1nextcheb(y, 3.30820231092092828324E-16, ref b0, ref b1, ref b2);
685                besselm1nextcheb(y, -1.88035477551078244854E-15, ref b0, ref b1, ref b2);
686                besselm1nextcheb(y, -3.81440307243700780478E-15, ref b0, ref b1, ref b2);
687                besselm1nextcheb(y, 1.04202769841288027642E-14, ref b0, ref b1, ref b2);
688                besselm1nextcheb(y, 4.27244001671195135429E-14, ref b0, ref b1, ref b2);
689                besselm1nextcheb(y, -2.10154184277266431302E-14, ref b0, ref b1, ref b2);
690                besselm1nextcheb(y, -4.08355111109219731823E-13, ref b0, ref b1, ref b2);
691                besselm1nextcheb(y, -7.19855177624590851209E-13, ref b0, ref b1, ref b2);
692                besselm1nextcheb(y, 2.03562854414708950722E-12, ref b0, ref b1, ref b2);
693                besselm1nextcheb(y, 1.41258074366137813316E-11, ref b0, ref b1, ref b2);
694                besselm1nextcheb(y, 3.25260358301548823856E-11, ref b0, ref b1, ref b2);
695                besselm1nextcheb(y, -1.89749581235054123450E-11, ref b0, ref b1, ref b2);
696                besselm1nextcheb(y, -5.58974346219658380687E-10, ref b0, ref b1, ref b2);
697                besselm1nextcheb(y, -3.83538038596423702205E-9, ref b0, ref b1, ref b2);
698                besselm1nextcheb(y, -2.63146884688951950684E-8, ref b0, ref b1, ref b2);
699                besselm1nextcheb(y, -2.51223623787020892529E-7, ref b0, ref b1, ref b2);
700                besselm1nextcheb(y, -3.88256480887769039346E-6, ref b0, ref b1, ref b2);
701                besselm1nextcheb(y, -1.10588938762623716291E-4, ref b0, ref b1, ref b2);
702                besselm1nextcheb(y, -9.76109749136146840777E-3, ref b0, ref b1, ref b2);
703                besselm1nextcheb(y, 7.78576235018280120474E-1, ref b0, ref b1, ref b2);
704                v = 0.5*(b0-b2);
705                z = v*Math.Exp(z)/Math.Sqrt(z);
706            }
707            if( (double)(x)<(double)(0) )
708            {
709                z = -z;
710            }
711            result = z;
712            return result;
713        }
714
715
716        /*************************************************************************
717        Modified Bessel function, second kind, order zero
718
719        Returns modified Bessel function of the second kind
720        of order zero of the argument.
721
722        The range is partitioned into the two intervals [0,8] and
723        (8, infinity).  Chebyshev polynomial expansions are employed
724        in each interval.
725
726        ACCURACY:
727
728        Tested at 2000 random points between 0 and 8.  Peak absolute
729        error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
730                             Relative error:
731        arithmetic   domain     # trials      peak         rms
732           IEEE      0, 30       30000       1.2e-15     1.6e-16
733
734        Cephes Math Library Release 2.8:  June, 2000
735        Copyright 1984, 1987, 2000 by Stephen L. Moshier
736        *************************************************************************/
737        public static double besselk0(double x)
738        {
739            double result = 0;
740            double y = 0;
741            double z = 0;
742            double v = 0;
743            double b0 = 0;
744            double b1 = 0;
745            double b2 = 0;
746
747            System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselK0: x<=0");
748            if( (double)(x)<=(double)(2) )
749            {
750                y = x*x-2.0;
751                besselmfirstcheb(1.37446543561352307156E-16, ref b0, ref b1, ref b2);
752                besselmnextcheb(y, 4.25981614279661018399E-14, ref b0, ref b1, ref b2);
753                besselmnextcheb(y, 1.03496952576338420167E-11, ref b0, ref b1, ref b2);
754                besselmnextcheb(y, 1.90451637722020886025E-9, ref b0, ref b1, ref b2);
755                besselmnextcheb(y, 2.53479107902614945675E-7, ref b0, ref b1, ref b2);
756                besselmnextcheb(y, 2.28621210311945178607E-5, ref b0, ref b1, ref b2);
757                besselmnextcheb(y, 1.26461541144692592338E-3, ref b0, ref b1, ref b2);
758                besselmnextcheb(y, 3.59799365153615016266E-2, ref b0, ref b1, ref b2);
759                besselmnextcheb(y, 3.44289899924628486886E-1, ref b0, ref b1, ref b2);
760                besselmnextcheb(y, -5.35327393233902768720E-1, ref b0, ref b1, ref b2);
761                v = 0.5*(b0-b2);
762                v = v-Math.Log(0.5*x)*besseli0(x);
763            }
764            else
765            {
766                z = 8.0/x-2.0;
767                besselmfirstcheb(5.30043377268626276149E-18, ref b0, ref b1, ref b2);
768                besselmnextcheb(z, -1.64758043015242134646E-17, ref b0, ref b1, ref b2);
769                besselmnextcheb(z, 5.21039150503902756861E-17, ref b0, ref b1, ref b2);
770                besselmnextcheb(z, -1.67823109680541210385E-16, ref b0, ref b1, ref b2);
771                besselmnextcheb(z, 5.51205597852431940784E-16, ref b0, ref b1, ref b2);
772                besselmnextcheb(z, -1.84859337734377901440E-15, ref b0, ref b1, ref b2);
773                besselmnextcheb(z, 6.34007647740507060557E-15, ref b0, ref b1, ref b2);
774                besselmnextcheb(z, -2.22751332699166985548E-14, ref b0, ref b1, ref b2);
775                besselmnextcheb(z, 8.03289077536357521100E-14, ref b0, ref b1, ref b2);
776                besselmnextcheb(z, -2.98009692317273043925E-13, ref b0, ref b1, ref b2);
777                besselmnextcheb(z, 1.14034058820847496303E-12, ref b0, ref b1, ref b2);
778                besselmnextcheb(z, -4.51459788337394416547E-12, ref b0, ref b1, ref b2);
779                besselmnextcheb(z, 1.85594911495471785253E-11, ref b0, ref b1, ref b2);
780                besselmnextcheb(z, -7.95748924447710747776E-11, ref b0, ref b1, ref b2);
781                besselmnextcheb(z, 3.57739728140030116597E-10, ref b0, ref b1, ref b2);
782                besselmnextcheb(z, -1.69753450938905987466E-9, ref b0, ref b1, ref b2);
783                besselmnextcheb(z, 8.57403401741422608519E-9, ref b0, ref b1, ref b2);
784                besselmnextcheb(z, -4.66048989768794782956E-8, ref b0, ref b1, ref b2);
785                besselmnextcheb(z, 2.76681363944501510342E-7, ref b0, ref b1, ref b2);
786                besselmnextcheb(z, -1.83175552271911948767E-6, ref b0, ref b1, ref b2);
787                besselmnextcheb(z, 1.39498137188764993662E-5, ref b0, ref b1, ref b2);
788                besselmnextcheb(z, -1.28495495816278026384E-4, ref b0, ref b1, ref b2);
789                besselmnextcheb(z, 1.56988388573005337491E-3, ref b0, ref b1, ref b2);
790                besselmnextcheb(z, -3.14481013119645005427E-2, ref b0, ref b1, ref b2);
791                besselmnextcheb(z, 2.44030308206595545468E0, ref b0, ref b1, ref b2);
792                v = 0.5*(b0-b2);
793                v = v*Math.Exp(-x)/Math.Sqrt(x);
794            }
795            result = v;
796            return result;
797        }
798
799
800        /*************************************************************************
801        Modified Bessel function, second kind, order one
802
803        Computes the modified Bessel function of the second kind
804        of order one of the argument.
805
806        The range is partitioned into the two intervals [0,2] and
807        (2, infinity).  Chebyshev polynomial expansions are employed
808        in each interval.
809
810        ACCURACY:
811
812                             Relative error:
813        arithmetic   domain     # trials      peak         rms
814           IEEE      0, 30       30000       1.2e-15     1.6e-16
815
816        Cephes Math Library Release 2.8:  June, 2000
817        Copyright 1984, 1987, 2000 by Stephen L. Moshier
818        *************************************************************************/
819        public static double besselk1(double x)
820        {
821            double result = 0;
822            double y = 0;
823            double z = 0;
824            double v = 0;
825            double b0 = 0;
826            double b1 = 0;
827            double b2 = 0;
828
829            z = 0.5*x;
830            System.Diagnostics.Debug.Assert((double)(z)>(double)(0), "Domain error in K1");
831            if( (double)(x)<=(double)(2) )
832            {
833                y = x*x-2.0;
834                besselm1firstcheb(-7.02386347938628759343E-18, ref b0, ref b1, ref b2);
835                besselm1nextcheb(y, -2.42744985051936593393E-15, ref b0, ref b1, ref b2);
836                besselm1nextcheb(y, -6.66690169419932900609E-13, ref b0, ref b1, ref b2);
837                besselm1nextcheb(y, -1.41148839263352776110E-10, ref b0, ref b1, ref b2);
838                besselm1nextcheb(y, -2.21338763073472585583E-8, ref b0, ref b1, ref b2);
839                besselm1nextcheb(y, -2.43340614156596823496E-6, ref b0, ref b1, ref b2);
840                besselm1nextcheb(y, -1.73028895751305206302E-4, ref b0, ref b1, ref b2);
841                besselm1nextcheb(y, -6.97572385963986435018E-3, ref b0, ref b1, ref b2);
842                besselm1nextcheb(y, -1.22611180822657148235E-1, ref b0, ref b1, ref b2);
843                besselm1nextcheb(y, -3.53155960776544875667E-1, ref b0, ref b1, ref b2);
844                besselm1nextcheb(y, 1.52530022733894777053E0, ref b0, ref b1, ref b2);
845                v = 0.5*(b0-b2);
846                result = Math.Log(z)*besseli1(x)+v/x;
847            }
848            else
849            {
850                y = 8.0/x-2.0;
851                besselm1firstcheb(-5.75674448366501715755E-18, ref b0, ref b1, ref b2);
852                besselm1nextcheb(y, 1.79405087314755922667E-17, ref b0, ref b1, ref b2);
853                besselm1nextcheb(y, -5.68946255844285935196E-17, ref b0, ref b1, ref b2);
854                besselm1nextcheb(y, 1.83809354436663880070E-16, ref b0, ref b1, ref b2);
855                besselm1nextcheb(y, -6.05704724837331885336E-16, ref b0, ref b1, ref b2);
856                besselm1nextcheb(y, 2.03870316562433424052E-15, ref b0, ref b1, ref b2);
857                besselm1nextcheb(y, -7.01983709041831346144E-15, ref b0, ref b1, ref b2);
858                besselm1nextcheb(y, 2.47715442448130437068E-14, ref b0, ref b1, ref b2);
859                besselm1nextcheb(y, -8.97670518232499435011E-14, ref b0, ref b1, ref b2);
860                besselm1nextcheb(y, 3.34841966607842919884E-13, ref b0, ref b1, ref b2);
861                besselm1nextcheb(y, -1.28917396095102890680E-12, ref b0, ref b1, ref b2);
862                besselm1nextcheb(y, 5.13963967348173025100E-12, ref b0, ref b1, ref b2);
863                besselm1nextcheb(y, -2.12996783842756842877E-11, ref b0, ref b1, ref b2);
864                besselm1nextcheb(y, 9.21831518760500529508E-11, ref b0, ref b1, ref b2);
865                besselm1nextcheb(y, -4.19035475934189648750E-10, ref b0, ref b1, ref b2);
866                besselm1nextcheb(y, 2.01504975519703286596E-9, ref b0, ref b1, ref b2);
867                besselm1nextcheb(y, -1.03457624656780970260E-8, ref b0, ref b1, ref b2);
868                besselm1nextcheb(y, 5.74108412545004946722E-8, ref b0, ref b1, ref b2);
869                besselm1nextcheb(y, -3.50196060308781257119E-7, ref b0, ref b1, ref b2);
870                besselm1nextcheb(y, 2.40648494783721712015E-6, ref b0, ref b1, ref b2);
871                besselm1nextcheb(y, -1.93619797416608296024E-5, ref b0, ref b1, ref b2);
872                besselm1nextcheb(y, 1.95215518471351631108E-4, ref b0, ref b1, ref b2);
873                besselm1nextcheb(y, -2.85781685962277938680E-3, ref b0, ref b1, ref b2);
874                besselm1nextcheb(y, 1.03923736576817238437E-1, ref b0, ref b1, ref b2);
875                besselm1nextcheb(y, 2.72062619048444266945E0, ref b0, ref b1, ref b2);
876                v = 0.5*(b0-b2);
877                result = Math.Exp(-x)*v/Math.Sqrt(x);
878            }
879            return result;
880        }
881
882
883        /*************************************************************************
884        Modified Bessel function, second kind, integer order
885
886        Returns modified Bessel function of the second kind
887        of order n of the argument.
888
889        The range is partitioned into the two intervals [0,9.55] and
890        (9.55, infinity).  An ascending power series is used in the
891        low range, and an asymptotic expansion in the high range.
892
893        ACCURACY:
894
895                             Relative error:
896        arithmetic   domain     # trials      peak         rms
897           IEEE      0,30        90000       1.8e-8      3.0e-10
898
899        Error is high only near the crossover point x = 9.55
900        between the two expansions used.
901
902        Cephes Math Library Release 2.8:  June, 2000
903        Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
904        *************************************************************************/
905        public static double besselkn(int nn,
906            double x)
907        {
908            double result = 0;
909            double k = 0;
910            double kf = 0;
911            double nk1f = 0;
912            double nkf = 0;
913            double zn = 0;
914            double t = 0;
915            double s = 0;
916            double z0 = 0;
917            double z = 0;
918            double ans = 0;
919            double fn = 0;
920            double pn = 0;
921            double pk = 0;
922            double zmn = 0;
923            double tlg = 0;
924            double tox = 0;
925            int i = 0;
926            int n = 0;
927            double eul = 0;
928
929            eul = 5.772156649015328606065e-1;
930            if( nn<0 )
931            {
932                n = -nn;
933            }
934            else
935            {
936                n = nn;
937            }
938            System.Diagnostics.Debug.Assert(n<=31, "Overflow in BesselKN");
939            System.Diagnostics.Debug.Assert((double)(x)>(double)(0), "Domain error in BesselKN");
940            if( (double)(x)<=(double)(9.55) )
941            {
942                ans = 0.0;
943                z0 = 0.25*x*x;
944                fn = 1.0;
945                pn = 0.0;
946                zmn = 1.0;
947                tox = 2.0/x;
948                if( n>0 )
949                {
950                    pn = -eul;
951                    k = 1.0;
952                    for(i=1; i<=n-1; i++)
953                    {
954                        pn = pn+1.0/k;
955                        k = k+1.0;
956                        fn = fn*k;
957                    }
958                    zmn = tox;
959                    if( n==1 )
960                    {
961                        ans = 1.0/x;
962                    }
963                    else
964                    {
965                        nk1f = fn/n;
966                        kf = 1.0;
967                        s = nk1f;
968                        z = -z0;
969                        zn = 1.0;
970                        for(i=1; i<=n-1; i++)
971                        {
972                            nk1f = nk1f/(n-i);
973                            kf = kf*i;
974                            zn = zn*z;
975                            t = nk1f*zn/kf;
976                            s = s+t;
977                            System.Diagnostics.Debug.Assert((double)(AP.Math.MaxRealNumber-Math.Abs(t))>(double)(Math.Abs(s)), "Overflow in BesselKN");
978                            System.Diagnostics.Debug.Assert(!((double)(tox)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/tox)<(double)(zmn)), "Overflow in BesselKN");
979                            zmn = zmn*tox;
980                        }
981                        s = s*0.5;
982                        t = Math.Abs(s);
983                        System.Diagnostics.Debug.Assert(!((double)(zmn)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/zmn)<(double)(t)), "Overflow in BesselKN");
984                        System.Diagnostics.Debug.Assert(!((double)(t)>(double)(1.0) & (double)(AP.Math.MaxRealNumber/t)<(double)(zmn)), "Overflow in BesselKN");
985                        ans = s*zmn;
986                    }
987                }
988                tlg = 2.0*Math.Log(0.5*x);
989                pk = -eul;
990                if( n==0 )
991                {
992                    pn = pk;
993                    t = 1.0;
994                }
995                else
996                {
997                    pn = pn+1.0/n;
998                    t = 1.0/fn;
999                }
1000                s = (pk+pn-tlg)*t;
1001                k = 1.0;
1002                do
1003                {
1004                    t = t*(z0/(k*(k+n)));
1005                    pk = pk+1.0/k;
1006                    pn = pn+1.0/(k+n);
1007                    s = s+(pk+pn-tlg)*t;
1008                    k = k+1.0;
1009                }
1010                while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) );
1011                s = 0.5*s/zmn;
1012                if( n%2!=0 )
1013                {
1014                    s = -s;
1015                }
1016                ans = ans+s;
1017                result = ans;
1018                return result;
1019            }
1020            if( (double)(x)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
1021            {
1022                result = 0;
1023                return result;
1024            }
1025            k = n;
1026            pn = 4.0*k*k;
1027            pk = 1.0;
1028            z0 = 8.0*x;
1029            fn = 1.0;
1030            t = 1.0;
1031            s = t;
1032            nkf = AP.Math.MaxRealNumber;
1033            i = 0;
1034            do
1035            {
1036                z = pn-pk*pk;
1037                t = t*z/(fn*z0);
1038                nk1f = Math.Abs(t);
1039                if( i>=n & (double)(nk1f)>(double)(nkf) )
1040                {
1041                    break;
1042                }
1043                nkf = nk1f;
1044                s = s+t;
1045                fn = fn+1.0;
1046                pk = pk+2.0;
1047                i = i+1;
1048            }
1049            while( (double)(Math.Abs(t/s))>(double)(AP.Math.MachineEpsilon) );
1050            result = Math.Exp(-x)*Math.Sqrt(Math.PI/(2.0*x))*s;
1051            return result;
1052        }
1053
1054
1055        /*************************************************************************
1056        Internal subroutine
1057
1058        Cephes Math Library Release 2.8:  June, 2000
1059        Copyright 1984, 1987, 2000 by Stephen L. Moshier
1060        *************************************************************************/
1061        private static void besselmfirstcheb(double c,
1062            ref double b0,
1063            ref double b1,
1064            ref double b2)
1065        {
1066            b0 = c;
1067            b1 = 0.0;
1068            b2 = 0.0;
1069        }
1070
1071
1072        /*************************************************************************
1073        Internal subroutine
1074
1075        Cephes Math Library Release 2.8:  June, 2000
1076        Copyright 1984, 1987, 2000 by Stephen L. Moshier
1077        *************************************************************************/
1078        private static void besselmnextcheb(double x,
1079            double c,
1080            ref double b0,
1081            ref double b1,
1082            ref double b2)
1083        {
1084            b2 = b1;
1085            b1 = b0;
1086            b0 = x*b1-b2+c;
1087        }
1088
1089
1090        /*************************************************************************
1091        Internal subroutine
1092
1093        Cephes Math Library Release 2.8:  June, 2000
1094        Copyright 1984, 1987, 2000 by Stephen L. Moshier
1095        *************************************************************************/
1096        private static void besselm1firstcheb(double c,
1097            ref double b0,
1098            ref double b1,
1099            ref double b2)
1100        {
1101            b0 = c;
1102            b1 = 0.0;
1103            b2 = 0.0;
1104        }
1105
1106
1107        /*************************************************************************
1108        Internal subroutine
1109
1110        Cephes Math Library Release 2.8:  June, 2000
1111        Copyright 1984, 1987, 2000 by Stephen L. Moshier
1112        *************************************************************************/
1113        private static void besselm1nextcheb(double x,
1114            double c,
1115            ref double b0,
1116            ref double b1,
1117            ref double b2)
1118        {
1119            b2 = b1;
1120            b1 = b0;
1121            b0 = x*b1-b2+c;
1122        }
1123
1124
1125        private static void besselasympt0(double x,
1126            ref double pzero,
1127            ref double qzero)
1128        {
1129            double xsq = 0;
1130            double p2 = 0;
1131            double q2 = 0;
1132            double p3 = 0;
1133            double q3 = 0;
1134
1135            xsq = 64.0/(x*x);
1136            p2 = 0.0;
1137            p2 = 2485.271928957404011288128951+xsq*p2;
1138            p2 = 153982.6532623911470917825993+xsq*p2;
1139            p2 = 2016135.283049983642487182349+xsq*p2;
1140            p2 = 8413041.456550439208464315611+xsq*p2;
1141            p2 = 12332384.76817638145232406055+xsq*p2;
1142            p2 = 5393485.083869438325262122897+xsq*p2;
1143            q2 = 1.0;
1144            q2 = 2615.700736920839685159081813+xsq*q2;
1145            q2 = 156001.7276940030940592769933+xsq*q2;
1146            q2 = 2025066.801570134013891035236+xsq*q2;
1147            q2 = 8426449.050629797331554404810+xsq*q2;
1148            q2 = 12338310.22786324960844856182+xsq*q2;
1149            q2 = 5393485.083869438325560444960+xsq*q2;
1150            p3 = -0.0;
1151            p3 = -4.887199395841261531199129300+xsq*p3;
1152            p3 = -226.2630641933704113967255053+xsq*p3;
1153            p3 = -2365.956170779108192723612816+xsq*p3;
1154            p3 = -8239.066313485606568803548860+xsq*p3;
1155            p3 = -10381.41698748464093880530341+xsq*p3;
1156            p3 = -3984.617357595222463506790588+xsq*p3;
1157            q3 = 1.0;
1158            q3 = 408.7714673983499223402830260+xsq*q3;
1159            q3 = 15704.89191515395519392882766+xsq*q3;
1160            q3 = 156021.3206679291652539287109+xsq*q3;
1161            q3 = 533291.3634216897168722255057+xsq*q3;
1162            q3 = 666745.4239319826986004038103+xsq*q3;
1163            q3 = 255015.5108860942382983170882+xsq*q3;
1164            pzero = p2/q2;
1165            qzero = 8*p3/q3/x;
1166        }
1167
1168
1169        private static void besselasympt1(double x,
1170            ref double pzero,
1171            ref double qzero)
1172        {
1173            double xsq = 0;
1174            double p2 = 0;
1175            double q2 = 0;
1176            double p3 = 0;
1177            double q3 = 0;
1178
1179            xsq = 64.0/(x*x);
1180            p2 = -1611.616644324610116477412898;
1181            p2 = -109824.0554345934672737413139+xsq*p2;
1182            p2 = -1523529.351181137383255105722+xsq*p2;
1183            p2 = -6603373.248364939109255245434+xsq*p2;
1184            p2 = -9942246.505077641195658377899+xsq*p2;
1185            p2 = -4435757.816794127857114720794+xsq*p2;
1186            q2 = 1.0;
1187            q2 = -1455.009440190496182453565068+xsq*q2;
1188            q2 = -107263.8599110382011903063867+xsq*q2;
1189            q2 = -1511809.506634160881644546358+xsq*q2;
1190            q2 = -6585339.479723087072826915069+xsq*q2;
1191            q2 = -9934124.389934585658967556309+xsq*q2;
1192            q2 = -4435757.816794127856828016962+xsq*q2;
1193            p3 = 35.26513384663603218592175580;
1194            p3 = 1706.375429020768002061283546+xsq*p3;
1195            p3 = 18494.26287322386679652009819+xsq*p3;
1196            p3 = 66178.83658127083517939992166+xsq*p3;
1197            p3 = 85145.16067533570196555001171+xsq*p3;
1198            p3 = 33220.91340985722351859704442+xsq*p3;
1199            q3 = 1.0;
1200            q3 = 863.8367769604990967475517183+xsq*q3;
1201            q3 = 37890.22974577220264142952256+xsq*q3;
1202            q3 = 400294.4358226697511708610813+xsq*q3;
1203            q3 = 1419460.669603720892855755253+xsq*q3;
1204            q3 = 1819458.042243997298924553839+xsq*q3;
1205            q3 = 708712.8194102874357377502472+xsq*q3;
1206            pzero = p2/q2;
1207            qzero = 8*p3/q3/x;
1208        }
1209    }
1210}
Note: See TracBrowser for help on using the repository browser.