source: sasview/src/sas/models/c_extension/cephes/j1.c @ 79492222

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 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.
  • Property mode set to 100644
File size: 11.8 KB
Line 
1/*                                                      j1.c
2 *
3 *      Bessel function of order one
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, j1();
10 *
11 * y = j1( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns Bessel function of order one of the argument.
18 *
19 * The domain is divided into the intervals [0, 8] and
20 * (8, infinity). In the first interval a 24 term Chebyshev
21 * expansion is used. In the second, the asymptotic
22 * trigonometric representation is employed using two
23 * rational functions of degree 5/5.
24 *
25 *
26 *
27 * ACCURACY:
28 *
29 *                      Absolute error:
30 * arithmetic   domain      # trials      peak         rms
31 *    DEC       0, 30       10000       4.0e-17     1.1e-17
32 *    IEEE      0, 30       30000       2.6e-16     1.1e-16
33 *
34 *
35 */
36/*                                                     y1.c
37 *
38 *      Bessel function of second kind of order one
39 *
40 *
41 *
42 * SYNOPSIS:
43 *
44 * double x, y, y1();
45 *
46 * y = y1( x );
47 *
48 *
49 *
50 * DESCRIPTION:
51 *
52 * Returns Bessel function of the second kind of order one
53 * of the argument.
54 *
55 * The domain is divided into the intervals [0, 8] and
56 * (8, infinity). In the first interval a 25 term Chebyshev
57 * expansion is used, and a call to j1() is required.
58 * In the second, the asymptotic trigonometric representation
59 * is employed using two rational functions of degree 5/5.
60 *
61 *
62 *
63 * ACCURACY:
64 *
65 *                      Absolute error:
66 * arithmetic   domain      # trials      peak         rms
67 *    DEC       0, 30       10000       8.6e-17     1.3e-17
68 *    IEEE      0, 30       30000       1.0e-15     1.3e-16
69 *
70 * (error criterion relative when |y1| > 1).
71 *
72 */
73
74
75/*
76Cephes Math Library Release 2.8:  June, 2000
77Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
78*/
79
80/*
81#define PIO4 .78539816339744830962
82#define THPIO4 2.35619449019234492885
83#define SQ2OPI .79788456080286535588
84*/
85
86#include "mconf.h"
87
88#ifdef UNK
89static double RP[4] = {
90-8.99971225705559398224E8,
91 4.52228297998194034323E11,
92-7.27494245221818276015E13,
93 3.68295732863852883286E15,
94};
95static double RQ[8] = {
96/* 1.00000000000000000000E0,*/
97 6.20836478118054335476E2,
98 2.56987256757748830383E5,
99 8.35146791431949253037E7,
100 2.21511595479792499675E10,
101 4.74914122079991414898E12,
102 7.84369607876235854894E14,
103 8.95222336184627338078E16,
104 5.32278620332680085395E18,
105};
106#endif
107#ifdef DEC
108static unsigned short RP[16] = {
1090147526,0110742,0063322,0077052,
1100051722,0112720,0065034,0061530,
1110153604,0052227,0033147,0105650,
1120055121,0055025,0032276,0022015,
113};
114static unsigned short RQ[32] = {
115/*0040200,0000000,0000000,0000000,*/
1160042433,0032610,0155604,0033473,
1170044572,0173320,0067270,0006616,
1180046637,0045246,0162225,0006606,
1190050645,0004773,0157577,0053004,
1200052612,0033734,0001667,0176501,
1210054462,0054121,0173147,0121367,
1220056237,0002777,0121451,0176007,
1230057623,0136253,0131601,0044710,
124};
125#endif
126#ifdef IBMPC
127static unsigned short RP[16] = {
1280x4fc5,0x4cda,0xd23c,0xc1ca,
1290x8c6b,0x0d43,0x52ba,0x425a,
1300xf175,0xe6cc,0x8a92,0xc2d0,
1310xc482,0xa697,0x2b42,0x432a,
132};
133static unsigned short RQ[32] = {
134/*0x0000,0x0000,0x0000,0x3ff0,*/
1350x86e7,0x1b70,0x66b1,0x4083,
1360x01b2,0x0dd7,0x5eda,0x410f,
1370xa1b1,0xdc92,0xe954,0x4193,
1380xeac1,0x7bef,0xa13f,0x4214,
1390xffa8,0x8076,0x46fb,0x4291,
1400xf45f,0x3ecc,0x4b0a,0x4306,
1410x3f81,0xf465,0xe0bf,0x4373,
1420x2939,0x7670,0x7795,0x43d2,
143};
144#endif
145#ifdef MIEEE
146static unsigned short RP[16] = {
1470xc1ca,0xd23c,0x4cda,0x4fc5,
1480x425a,0x52ba,0x0d43,0x8c6b,
1490xc2d0,0x8a92,0xe6cc,0xf175,
1500x432a,0x2b42,0xa697,0xc482,
151};
152static unsigned short RQ[32] = {
153/*0x3ff0,0x0000,0x0000,0x0000,*/
1540x4083,0x66b1,0x1b70,0x86e7,
1550x410f,0x5eda,0x0dd7,0x01b2,
1560x4193,0xe954,0xdc92,0xa1b1,
1570x4214,0xa13f,0x7bef,0xeac1,
1580x4291,0x46fb,0x8076,0xffa8,
1590x4306,0x4b0a,0x3ecc,0xf45f,
1600x4373,0xe0bf,0xf465,0x3f81,
1610x43d2,0x7795,0x7670,0x2939,
162};
163#endif
164
165#ifdef UNK
166static double PP[7] = {
167 7.62125616208173112003E-4,
168 7.31397056940917570436E-2,
169 1.12719608129684925192E0,
170 5.11207951146807644818E0,
171 8.42404590141772420927E0,
172 5.21451598682361504063E0,
173 1.00000000000000000254E0,
174};
175static double PQ[7] = {
176 5.71323128072548699714E-4,
177 6.88455908754495404082E-2,
178 1.10514232634061696926E0,
179 5.07386386128601488557E0,
180 8.39985554327604159757E0,
181 5.20982848682361821619E0,
182 9.99999999999999997461E-1,
183};
184#endif
185#ifdef DEC
186static unsigned short PP[28] = {
1870035507,0144542,0061543,0024326,
1880037225,0145105,0017766,0022661,
1890040220,0043766,0010254,0133255,
1900040643,0113047,0142611,0151521,
1910041006,0144344,0055351,0074261,
1920040646,0156520,0120574,0006416,
1930040200,0000000,0000000,0000000,
194};
195static unsigned short PQ[28] = {
1960035425,0142330,0115041,0165514,
1970037214,0177352,0145105,0052026,
1980040215,0072515,0141207,0073255,
1990040642,0056427,0137222,0106405,
2000041006,0062716,0166427,0165450,
2010040646,0133352,0035425,0123304,
2020040200,0000000,0000000,0000000,
203};
204#endif
205#ifdef IBMPC
206static unsigned short PP[28] = {
2070x651b,0x4c6c,0xf92c,0x3f48,
2080xc4b6,0xa3fe,0xb948,0x3fb2,
2090x96d6,0xc215,0x08fe,0x3ff2,
2100x3a6a,0xf8b1,0x72c4,0x4014,
2110x2f16,0x8b5d,0xd91c,0x4020,
2120x81a2,0x142f,0xdbaa,0x4014,
2130x0000,0x0000,0x0000,0x3ff0,
214};
215static unsigned short PQ[28] = {
2160x3d69,0x1344,0xb89b,0x3f42,
2170xaa83,0x5948,0x9fdd,0x3fb1,
2180xeed6,0xb850,0xaea9,0x3ff1,
2190x51a1,0xf7d2,0x4ba2,0x4014,
2200xfd65,0xdda2,0xccb9,0x4020,
2210xb4d9,0x4762,0xd6dd,0x4014,
2220x0000,0x0000,0x0000,0x3ff0,
223};
224#endif
225#ifdef MIEEE
226static unsigned short PP[28] = {
2270x3f48,0xf92c,0x4c6c,0x651b,
2280x3fb2,0xb948,0xa3fe,0xc4b6,
2290x3ff2,0x08fe,0xc215,0x96d6,
2300x4014,0x72c4,0xf8b1,0x3a6a,
2310x4020,0xd91c,0x8b5d,0x2f16,
2320x4014,0xdbaa,0x142f,0x81a2,
2330x3ff0,0x0000,0x0000,0x0000,
234};
235static unsigned short PQ[28] = {
2360x3f42,0xb89b,0x1344,0x3d69,
2370x3fb1,0x9fdd,0x5948,0xaa83,
2380x3ff1,0xaea9,0xb850,0xeed6,
2390x4014,0x4ba2,0xf7d2,0x51a1,
2400x4020,0xccb9,0xdda2,0xfd65,
2410x4014,0xd6dd,0x4762,0xb4d9,
2420x3ff0,0x0000,0x0000,0x0000,
243};
244#endif
245
246#ifdef UNK
247static double QP[8] = {
248 5.10862594750176621635E-2,
249 4.98213872951233449420E0,
250 7.58238284132545283818E1,
251 3.66779609360150777800E2,
252 7.10856304998926107277E2,
253 5.97489612400613639965E2,
254 2.11688757100572135698E2,
255 2.52070205858023719784E1,
256};
257static double QQ[7] = {
258/* 1.00000000000000000000E0,*/
259 7.42373277035675149943E1,
260 1.05644886038262816351E3,
261 4.98641058337653607651E3,
262 9.56231892404756170795E3,
263 7.99704160447350683650E3,
264 2.82619278517639096600E3,
265 3.36093607810698293419E2,
266};
267#endif
268#ifdef DEC
269static unsigned short QP[32] = {
2700037121,0037723,0055605,0151004,
2710040637,0066656,0031554,0077264,
2720041627,0122714,0153170,0161466,
2730042267,0061712,0036520,0140145,
2740042461,0133315,0131573,0071176,
2750042425,0057525,0147500,0013201,
2760042123,0130122,0061245,0154131,
2770041311,0123772,0064254,0172650,
278};
279static unsigned short QQ[28] = {
280/*0040200,0000000,0000000,0000000,*/
2810041624,0074603,0002112,0101670,
2820042604,0007135,0010162,0175565,
2830043233,0151510,0157757,0172010,
2840043425,0064506,0112006,0104276,
2850043371,0164125,0032271,0164242,
2860043060,0121425,0122750,0136013,
2870042250,0005773,0053472,0146267,
288};
289#endif
290#ifdef IBMPC
291static unsigned short QP[32] = {
2920xba40,0x6b70,0x27fa,0x3faa,
2930x8fd6,0xc66d,0xedb5,0x4013,
2940x1c67,0x9acf,0xf4b9,0x4052,
2950x180d,0x47aa,0xec79,0x4076,
2960x6e50,0xb66f,0x36d9,0x4086,
2970x02d0,0xb9e8,0xabea,0x4082,
2980xbb0b,0x4c54,0x760a,0x406a,
2990x9eb5,0x4d15,0x34ff,0x4039,
300};
301static unsigned short QQ[28] = {
302/*0x0000,0x0000,0x0000,0x3ff0,*/
3030x5077,0x6089,0x8f30,0x4052,
3040x5f6f,0xa20e,0x81cb,0x4090,
3050xfe81,0x1bfd,0x7a69,0x40b3,
3060xd118,0xd280,0xad28,0x40c2,
3070x3d14,0xa697,0x3d0a,0x40bf,
3080x1781,0xb4bd,0x1462,0x40a6,
3090x5997,0x6ae7,0x017f,0x4075,
310};
311#endif
312#ifdef MIEEE
313static unsigned short QP[32] = {
3140x3faa,0x27fa,0x6b70,0xba40,
3150x4013,0xedb5,0xc66d,0x8fd6,
3160x4052,0xf4b9,0x9acf,0x1c67,
3170x4076,0xec79,0x47aa,0x180d,
3180x4086,0x36d9,0xb66f,0x6e50,
3190x4082,0xabea,0xb9e8,0x02d0,
3200x406a,0x760a,0x4c54,0xbb0b,
3210x4039,0x34ff,0x4d15,0x9eb5,
322};
323static unsigned short QQ[28] = {
324/*0x3ff0,0x0000,0x0000,0x0000,*/
3250x4052,0x8f30,0x6089,0x5077,
3260x4090,0x81cb,0xa20e,0x5f6f,
3270x40b3,0x7a69,0x1bfd,0xfe81,
3280x40c2,0xad28,0xd280,0xd118,
3290x40bf,0x3d0a,0xa697,0x3d14,
3300x40a6,0x1462,0xb4bd,0x1781,
3310x4075,0x017f,0x6ae7,0x5997,
332};
333#endif
334
335#ifdef UNK
336static double YP[6] = {
337 1.26320474790178026440E9,
338-6.47355876379160291031E11,
339 1.14509511541823727583E14,
340-8.12770255501325109621E15,
341 2.02439475713594898196E17,
342-7.78877196265950026825E17,
343};
344static double YQ[8] = {
345/* 1.00000000000000000000E0,*/
346 5.94301592346128195359E2,
347 2.35564092943068577943E5,
348 7.34811944459721705660E7,
349 1.87601316108706159478E10,
350 3.88231277496238566008E12,
351 6.20557727146953693363E14,
352 6.87141087355300489866E16,
353 3.97270608116560655612E18,
354};
355#endif
356#ifdef DEC
357static unsigned short YP[24] = {
3580047626,0112763,0013715,0133045,
3590152026,0134552,0142033,0024411,
3600053720,0045245,0102210,0077565,
3610155347,0000321,0136415,0102031,
3620056463,0146550,0055633,0032605,
3630157054,0171012,0167361,0054265,
364};
365static unsigned short YQ[32] = {
366/*0040200,0000000,0000000,0000000,*/
3670042424,0111515,0044773,0153014,
3680044546,0005405,0171307,0075774,
3690046614,0023575,0047105,0063556,
3700050613,0143034,0101533,0156026,
3710052541,0175367,0166514,0114257,
3720054415,0014466,0134350,0171154,
3730056164,0017436,0025075,0022101,
3740057534,0103614,0103663,0121772,
375};
376#endif
377#ifdef IBMPC
378static unsigned short YP[24] = {
3790xb6c5,0x62f9,0xd2be,0x41d2,
3800x6521,0x5883,0xd72d,0xc262,
3810x0fef,0xb091,0x0954,0x42da,
3820xb083,0x37a1,0xe01a,0xc33c,
3830x66b1,0x0b73,0x79ad,0x4386,
3840x2b17,0x5dde,0x9e41,0xc3a5,
385};
386static unsigned short YQ[32] = {
387/*0x0000,0x0000,0x0000,0x3ff0,*/
3880x7ac2,0xa93f,0x9269,0x4082,
3890xef7f,0xbe58,0xc160,0x410c,
3900xacee,0xa9c8,0x84ef,0x4191,
3910x7b83,0x906b,0x78c3,0x4211,
3920x9316,0xfda9,0x3f5e,0x428c,
3930x1e4e,0xd71d,0xa326,0x4301,
3940xa488,0xc547,0x83e3,0x436e,
3950x747f,0x90f6,0x90f1,0x43cb,
396};
397#endif
398#ifdef MIEEE
399static unsigned short YP[24] = {
4000x41d2,0xd2be,0x62f9,0xb6c5,
4010xc262,0xd72d,0x5883,0x6521,
4020x42da,0x0954,0xb091,0x0fef,
4030xc33c,0xe01a,0x37a1,0xb083,
4040x4386,0x79ad,0x0b73,0x66b1,
4050xc3a5,0x9e41,0x5dde,0x2b17,
406};
407static unsigned short YQ[32] = {
408/*0x3ff0,0x0000,0x0000,0x0000,*/
4090x4082,0x9269,0xa93f,0x7ac2,
4100x410c,0xc160,0xbe58,0xef7f,
4110x4191,0x84ef,0xa9c8,0xacee,
4120x4211,0x78c3,0x906b,0x7b83,
4130x428c,0x3f5e,0xfda9,0x9316,
4140x4301,0xa326,0xd71d,0x1e4e,
4150x436e,0x83e3,0xc547,0xa488,
4160x43cb,0x90f1,0x90f6,0x747f,
417};
418#endif
419
420
421#ifdef UNK
422static double Z1 = 1.46819706421238932572E1;
423static double Z2 = 4.92184563216946036703E1;
424#endif
425
426#ifdef DEC
427static unsigned short DZ1[] = {0041152,0164532,0006114,0010540};
428static unsigned short DZ2[] = {0041504,0157663,0001625,0020621};
429#define Z1 (*(double *)DZ1)
430#define Z2 (*(double *)DZ2)
431#endif
432
433#ifdef IBMPC
434static unsigned short DZ1[] = {0x822c,0x4189,0x5d2b,0x402d};
435static unsigned short DZ2[] = {0xa432,0x6072,0x9bf6,0x4048};
436#define Z1 (*(double *)DZ1)
437#define Z2 (*(double *)DZ2)
438#endif
439
440#ifdef MIEEE
441static unsigned short DZ1[] = {0x402d,0x5d2b,0x4189,0x822c};
442static unsigned short DZ2[] = {0x4048,0x9bf6,0x6072,0xa432};
443#define Z1 (*(double *)DZ1)
444#define Z2 (*(double *)DZ2)
445#endif
446
447#ifdef ANSIPROT
448extern double polevl ( double, void *, int );
449extern double p1evl ( double, void *, int );
450extern double log ( double );
451extern double sin ( double );
452extern double cos ( double );
453extern double sqrt ( double );
454double j1 ( double );
455#else
456double polevl(), p1evl(), log(), sin(), cos(), sqrt();
457double j1();
458#endif
459extern double TWOOPI, THPIO4, SQ2OPI;
460
461double j1(x)
462double x;
463{
464double w, z, p, q, xn;
465
466w = x;
467if( x < 0 )
468        w = -x;
469
470if( w <= 5.0 )
471        {
472        z = x * x;     
473        w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
474        w = w * x * (z - Z1) * (z - Z2);
475        return( w );
476        }
477
478w = 5.0/x;
479z = w * w;
480p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
481q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
482xn = x - THPIO4;
483p = p * cos(xn) - w * q * sin(xn);
484return( p * SQ2OPI / sqrt(x) );
485}
486
487
488extern double MAXNUM;
489
490double y1(x)
491double x;
492{
493double w, z, p, q, xn;
494
495if( x <= 5.0 )
496        {
497        if( x <= 0.0 )
498                {
499                mtherr( "y1", DOMAIN );
500                return( -MAXNUM );
501                }
502        z = x * x;
503        w = x * (polevl( z, YP, 5 ) / p1evl( z, YQ, 8 ));
504        w += TWOOPI * ( j1(x) * log(x)  -  1.0/x );
505        return( w );
506        }
507
508w = 5.0/x;
509z = w * w;
510p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
511q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
512xn = x - THPIO4;
513p = p * sin(xn) + w * q * cos(xn);
514return( p * SQ2OPI / sqrt(x) );
515}
Note: See TracBrowser for help on using the repository browser.