source: sasview/src/sas/models/c_extension/cephes/k1.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: 7.1 KB
Line 
1/*                                                      k1.c
2 *
3 *      Modified Bessel function, third kind, order one
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, k1();
10 *
11 * y = k1( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Computes the modified Bessel function of the third kind
18 * of order one of the argument.
19 *
20 * The range is partitioned into the two intervals [0,2] and
21 * (2, infinity).  Chebyshev polynomial expansions are employed
22 * in each interval.
23 *
24 *
25 *
26 * ACCURACY:
27 *
28 *                      Relative error:
29 * arithmetic   domain     # trials      peak         rms
30 *    DEC       0, 30        3300       8.9e-17     2.2e-17
31 *    IEEE      0, 30       30000       1.2e-15     1.6e-16
32 *
33 * ERROR MESSAGES:
34 *
35 *   message         condition      value returned
36 * k1 domain          x <= 0          MAXNUM
37 *
38 */
39/*                                                     k1e.c
40 *
41 *      Modified Bessel function, third kind, order one,
42 *      exponentially scaled
43 *
44 *
45 *
46 * SYNOPSIS:
47 *
48 * double x, y, k1e();
49 *
50 * y = k1e( x );
51 *
52 *
53 *
54 * DESCRIPTION:
55 *
56 * Returns exponentially scaled modified Bessel function
57 * of the third kind of order one of the argument:
58 *
59 *      k1e(x) = exp(x) * k1(x).
60 *
61 *
62 *
63 * ACCURACY:
64 *
65 *                      Relative error:
66 * arithmetic   domain     # trials      peak         rms
67 *    IEEE      0, 30       30000       7.8e-16     1.2e-16
68 * See k1().
69 *
70 */
71
72/*
73Cephes Math Library Release 2.8:  June, 2000
74Copyright 1984, 1987, 2000 by Stephen L. Moshier
75*/
76
77#include "mconf.h"
78
79/* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
80 * in the interval [0,2].
81 *
82 * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
83 */
84
85#ifdef UNK
86static double A[] =
87{
88-7.02386347938628759343E-18,
89-2.42744985051936593393E-15,
90-6.66690169419932900609E-13,
91-1.41148839263352776110E-10,
92-2.21338763073472585583E-8,
93-2.43340614156596823496E-6,
94-1.73028895751305206302E-4,
95-6.97572385963986435018E-3,
96-1.22611180822657148235E-1,
97-3.53155960776544875667E-1,
98 1.52530022733894777053E0
99};
100#endif
101
102#ifdef DEC
103static unsigned short A[] = {
1040122001,0110501,0164746,0151255,
1050124056,0165213,0150034,0147377,
1060126073,0124026,0167207,0001044,
1070130033,0030735,0141061,0033116,
1080131676,0020350,0121341,0107175,
1090133443,0046631,0062031,0070716,
1100135065,0067427,0026435,0164022,
1110136344,0112234,0165752,0006222,
1120137373,0015622,0017016,0155636,
1130137664,0150333,0125730,0067240,
1140040303,0036411,0130200,0043120
115};
116#endif
117
118#ifdef IBMPC
119static unsigned short A[] = {
1200xda56,0x3d3c,0x3228,0xbc60,
1210x99e0,0x7a03,0xdd51,0xbce5,
1220xe045,0xddd0,0x7502,0xbd67,
1230x26ca,0xb846,0x663b,0xbde3,
1240x31d0,0x145c,0xc41d,0xbe57,
1250x2e3a,0x2c83,0x69b3,0xbec4,
1260xbd02,0xe5a3,0xade2,0xbf26,
1270x4192,0x9d7d,0x9293,0xbf7c,
1280xdb74,0x43c1,0x6372,0xbfbf,
1290x0dd4,0x757b,0x9a1b,0xbfd6,
1300x08ca,0x3610,0x67a1,0x3ff8
131};
132#endif
133
134#ifdef MIEEE
135static unsigned short A[] = {
1360xbc60,0x3228,0x3d3c,0xda56,
1370xbce5,0xdd51,0x7a03,0x99e0,
1380xbd67,0x7502,0xddd0,0xe045,
1390xbde3,0x663b,0xb846,0x26ca,
1400xbe57,0xc41d,0x145c,0x31d0,
1410xbec4,0x69b3,0x2c83,0x2e3a,
1420xbf26,0xade2,0xe5a3,0xbd02,
1430xbf7c,0x9293,0x9d7d,0x4192,
1440xbfbf,0x6372,0x43c1,0xdb74,
1450xbfd6,0x9a1b,0x757b,0x0dd4,
1460x3ff8,0x67a1,0x3610,0x08ca
147};
148#endif
149
150
151
152/* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
153 * in the interval [2,infinity].
154 *
155 * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
156 */
157
158#ifdef UNK
159static double B[] =
160{
161-5.75674448366501715755E-18,
162 1.79405087314755922667E-17,
163-5.68946255844285935196E-17,
164 1.83809354436663880070E-16,
165-6.05704724837331885336E-16,
166 2.03870316562433424052E-15,
167-7.01983709041831346144E-15,
168 2.47715442448130437068E-14,
169-8.97670518232499435011E-14,
170 3.34841966607842919884E-13,
171-1.28917396095102890680E-12,
172 5.13963967348173025100E-12,
173-2.12996783842756842877E-11,
174 9.21831518760500529508E-11,
175-4.19035475934189648750E-10,
176 2.01504975519703286596E-9,
177-1.03457624656780970260E-8,
178 5.74108412545004946722E-8,
179-3.50196060308781257119E-7,
180 2.40648494783721712015E-6,
181-1.93619797416608296024E-5,
182 1.95215518471351631108E-4,
183-2.85781685962277938680E-3,
184 1.03923736576817238437E-1,
185 2.72062619048444266945E0
186};
187#endif
188
189#ifdef DEC
190static unsigned short B[] = {
1910121724,0061352,0013041,0150076,
1920022245,0074324,0016172,0173232,
1930122603,0030250,0135670,0165221,
1940023123,0165362,0023561,0060124,
1950123456,0112436,0141654,0073623,
1960024022,0163557,0077564,0006753,
1970124374,0165221,0131014,0026524,
1980024737,0017512,0144250,0175451,
1990125312,0021456,0123136,0076633,
2000025674,0077720,0020125,0102607,
2010126265,0067543,0007744,0043701,
2020026664,0152702,0033002,0074202,
2030127273,0055234,0120016,0071733,
2040027712,0133200,0042441,0075515,
2050130346,0057000,0015456,0074470,
2060031012,0074441,0051636,0111155,
2070131461,0136444,0177417,0002101,
2080032166,0111743,0032176,0021410,
2090132674,0001224,0076555,0027060,
2100033441,0077430,0135226,0106663,
2110134242,0065610,0167155,0113447,
2120035114,0131304,0043664,0102163,
2130136073,0045065,0171465,0122123,
2140037324,0152767,0147401,0017732,
2150040456,0017275,0050061,0062120,
216};
217#endif
218
219#ifdef IBMPC
220static unsigned short B[] = {
2210x3a08,0x42c4,0x8c5d,0xbc5a,
2220x5ed3,0x838f,0xaf1a,0x3c74,
2230x1d52,0x1777,0x6615,0xbc90,
2240x2c0b,0x44ee,0x7d5e,0x3caa,
2250x8ef2,0xd875,0xd2a3,0xbcc5,
2260x81bd,0xefee,0x5ced,0x3ce2,
2270x85ab,0x3641,0x9d52,0xbcff,
2280x1f65,0x5915,0xe3e9,0x3d1b,
2290xcfb3,0xd4cb,0x4465,0xbd39,
2300xb0b1,0x040a,0x8ffa,0x3d57,
2310x88f8,0x61fc,0xadec,0xbd76,
2320x4f10,0x46c0,0x9ab8,0x3d96,
2330xce7b,0x9401,0x6b53,0xbdb7,
2340x2f6a,0x08a4,0x56d0,0x3dd9,
2350xcf27,0x0365,0xcbc0,0xbdfc,
2360xd24e,0x2a73,0x4f24,0x3e21,
2370xe088,0x9fe1,0x37a4,0xbe46,
2380xc461,0x668f,0xd27c,0x3e6e,
2390xa5c6,0x8fad,0x8052,0xbe97,
2400xd1b6,0x1752,0x2fe3,0x3ec4,
2410xb2e5,0x1dcd,0x4d71,0xbef4,
2420x908e,0x88f6,0x9658,0x3f29,
2430xb48a,0xbe66,0x6946,0xbf67,
2440x23fb,0xf9e0,0x9abe,0x3fba,
2450x2c8a,0xaa06,0xc3d7,0x4005
246};
247#endif
248
249#ifdef MIEEE
250static unsigned short B[] = {
2510xbc5a,0x8c5d,0x42c4,0x3a08,
2520x3c74,0xaf1a,0x838f,0x5ed3,
2530xbc90,0x6615,0x1777,0x1d52,
2540x3caa,0x7d5e,0x44ee,0x2c0b,
2550xbcc5,0xd2a3,0xd875,0x8ef2,
2560x3ce2,0x5ced,0xefee,0x81bd,
2570xbcff,0x9d52,0x3641,0x85ab,
2580x3d1b,0xe3e9,0x5915,0x1f65,
2590xbd39,0x4465,0xd4cb,0xcfb3,
2600x3d57,0x8ffa,0x040a,0xb0b1,
2610xbd76,0xadec,0x61fc,0x88f8,
2620x3d96,0x9ab8,0x46c0,0x4f10,
2630xbdb7,0x6b53,0x9401,0xce7b,
2640x3dd9,0x56d0,0x08a4,0x2f6a,
2650xbdfc,0xcbc0,0x0365,0xcf27,
2660x3e21,0x4f24,0x2a73,0xd24e,
2670xbe46,0x37a4,0x9fe1,0xe088,
2680x3e6e,0xd27c,0x668f,0xc461,
2690xbe97,0x8052,0x8fad,0xa5c6,
2700x3ec4,0x2fe3,0x1752,0xd1b6,
2710xbef4,0x4d71,0x1dcd,0xb2e5,
2720x3f29,0x9658,0x88f6,0x908e,
2730xbf67,0x6946,0xbe66,0xb48a,
2740x3fba,0x9abe,0xf9e0,0x23fb,
2750x4005,0xc3d7,0xaa06,0x2c8a
276};
277#endif
278
279#ifdef ANSIPROT
280extern double chbevl ( double, void *, int );
281extern double exp ( double );
282extern double i1 ( double );
283extern double log ( double );
284extern double sqrt ( double );
285#else
286double chbevl(), exp(), i1(), log(), sqrt();
287#endif
288extern double PI;
289extern double MINLOG, MAXNUM;
290
291double k1(x)
292double x;
293{
294double y, z;
295
296z = 0.5 * x;
297if( z <= 0.0 )
298        {
299        mtherr( "k1", DOMAIN );
300        return( MAXNUM );
301        }
302
303if( x <= 2.0 )
304        {
305        y = x * x - 2.0;
306        y =  log(z) * i1(x)  +  chbevl( y, A, 11 ) / x;
307        return( y );
308        }
309
310return(  exp(-x) * chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x) );
311}
312
313
314
315
316double k1e( x )
317double x;
318{
319double y;
320
321if( x <= 0.0 )
322        {
323        mtherr( "k1e", DOMAIN );
324        return( MAXNUM );
325        }
326
327if( x <= 2.0 )
328        {
329        y = x * x - 2.0;
330        y =  log( 0.5 * x ) * i1(x)  +  chbevl( y, A, 11 ) / x;
331        return( y * exp(x) );
332        }
333
334return(  chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x) );
335}
Note: See TracBrowser for help on using the repository browser.