source: sasview/src/sas/models/c_extension/cephes/kn.c @ 7bcbeb7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 7bcbeb7 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 4.2 KB
Line 
1/*                                                      kn.c
2 *
3 *      Modified Bessel function, third kind, integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, kn();
10 * int n;
11 *
12 * y = kn( n, x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns modified Bessel function of the third kind
19 * of order n of the argument.
20 *
21 * The range is partitioned into the two intervals [0,9.55] and
22 * (9.55, infinity).  An ascending power series is used in the
23 * low range, and an asymptotic expansion in the high range.
24 *
25 *
26 *
27 * ACCURACY:
28 *
29 *                      Relative error:
30 * arithmetic   domain     # trials      peak         rms
31 *    DEC       0,30         3000       1.3e-9      5.8e-11
32 *    IEEE      0,30        90000       1.8e-8      3.0e-10
33 *
34 *  Error is high only near the crossover point x = 9.55
35 * between the two expansions used.
36 */
37
38
39/*
40Cephes Math Library Release 2.8:  June, 2000
41Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
42*/
43
44
45/*
46Algorithm for Kn.
47                       n-1
48                   -n   -  (n-k-1)!    2   k
49K (x)  =  0.5 (x/2)     >  -------- (-x /4)
50 n                      -     k!
51                       k=0
52
53                    inf.                                   2   k
54       n         n   -                                   (x /4)
55 + (-1)  0.5(x/2)    >  {p(k+1) + p(n+k+1) - 2log(x/2)} ---------
56                     -                                  k! (n+k)!
57                    k=0
58
59where  p(m) is the psi function: p(1) = -EUL and
60
61                      m-1
62                       -
63      p(m)  =  -EUL +  >  1/k
64                       -
65                      k=1
66
67For large x,
68                                         2        2     2
69                                      u-1     (u-1 )(u-3 )
70K (z)  =  sqrt(pi/2z) exp(-z) { 1 + ------- + ------------ + ...}
71 v                                        1            2
72                                    1! (8z)     2! (8z)
73asymptotically, where
74
75           2
76    u = 4 v .
77
78*/
79
80#include "mconf.h"
81
82#define EUL 5.772156649015328606065e-1
83#define MAXFAC 31
84#ifdef ANSIPROT
85extern double fabs ( double );
86extern double exp ( double );
87extern double log ( double );
88extern double sqrt ( double );
89#else
90double fabs(), exp(), log(), sqrt();
91#endif
92extern double MACHEP, MAXNUM, MAXLOG, PI;
93
94double kn( nn, x )
95int nn;
96double x;
97{
98double k, kf, nk1f, nkf, zn, t, s, z0, z;
99double ans, fn, pn, pk, zmn, tlg, tox;
100int i, n;
101
102if( nn < 0 )
103        n = -nn;
104else
105        n = nn;
106
107if( n > MAXFAC )
108        {
109overf:
110        mtherr( "kn", OVERFLOW );
111        return( MAXNUM );
112        }
113
114if( x <= 0.0 )
115        {
116        if( x < 0.0 )
117                mtherr( "kn", DOMAIN );
118        else
119                mtherr( "kn", SING );
120        return( MAXNUM );
121        }
122
123
124if( x > 9.55 )
125        goto asymp;
126
127ans = 0.0;
128z0 = 0.25 * x * x;
129fn = 1.0;
130pn = 0.0;
131zmn = 1.0;
132tox = 2.0/x;
133
134if( n > 0 )
135        {
136        /* compute factorial of n and psi(n) */
137        pn = -EUL;
138        k = 1.0;
139        for( i=1; i<n; i++ )
140                {
141                pn += 1.0/k;
142                k += 1.0;
143                fn *= k;
144                }
145
146        zmn = tox;
147
148        if( n == 1 )
149                {
150                ans = 1.0/x;
151                }
152        else
153                {
154                nk1f = fn/n;
155                kf = 1.0;
156                s = nk1f;
157                z = -z0;
158                zn = 1.0;
159                for( i=1; i<n; i++ )
160                        {
161                        nk1f = nk1f/(n-i);
162                        kf = kf * i;
163                        zn *= z;
164                        t = nk1f * zn / kf;
165                        s += t;   
166                        if( (MAXNUM - fabs(t)) < fabs(s) )
167                                goto overf;
168                        if( (tox > 1.0) && ((MAXNUM/tox) < zmn) )
169                                goto overf;
170                        zmn *= tox;
171                        }
172                s *= 0.5;
173                t = fabs(s);
174                if( (zmn > 1.0) && ((MAXNUM/zmn) < t) )
175                        goto overf;
176                if( (t > 1.0) && ((MAXNUM/t) < zmn) )
177                        goto overf;
178                ans = s * zmn;
179                }
180        }
181
182
183tlg = 2.0 * log( 0.5 * x );
184pk = -EUL;
185if( n == 0 )
186        {
187        pn = pk;
188        t = 1.0;
189        }
190else
191        {
192        pn = pn + 1.0/n;
193        t = 1.0/fn;
194        }
195s = (pk+pn-tlg)*t;
196k = 1.0;
197do
198        {
199        t *= z0 / (k * (k+n));
200        pk += 1.0/k;
201        pn += 1.0/(k+n);
202        s += (pk+pn-tlg)*t;
203        k += 1.0;
204        }
205while( fabs(t/s) > MACHEP );
206
207s = 0.5 * s / zmn;
208if( n & 1 )
209        s = -s;
210ans += s;
211
212return(ans);
213
214
215
216/* Asymptotic expansion for Kn(x) */
217/* Converges to 1.4e-17 for x > 18.4 */
218
219asymp:
220
221if( x > MAXLOG )
222        {
223        mtherr( "kn", UNDERFLOW );
224        return(0.0);
225        }
226k = n;
227pn = 4.0 * k * k;
228pk = 1.0;
229z0 = 8.0 * x;
230fn = 1.0;
231t = 1.0;
232s = t;
233nkf = MAXNUM;
234i = 0;
235do
236        {
237        z = pn - pk * pk;
238        t = t * z /(fn * z0);
239        nk1f = fabs(t);
240        if( (i >= n) && (nk1f > nkf) )
241                {
242                goto adone;
243                }
244        nkf = nk1f;
245        s += t;
246        fn += 1.0;
247        pk += 2.0;
248        i += 1;
249        }
250while( fabs(t/s) > MACHEP );
251
252adone:
253ans = exp(-x) * sqrt( PI/(2.0*x) ) * s;
254return(ans);
255}
Note: See TracBrowser for help on using the repository browser.