source: sasmodels/sasmodels/models/lib/j1d.c @ 3936ad3

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3936ad3 was 3936ad3, checked in by wojciech, 8 years ago

Bessel functions from double-precison cephes has been added

  • Property mode set to 100644
File size: 4.0 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/*
75Cephes Math Library Release 2.8:  June, 2000
76Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
77*/
78double j1(double );
79
80double j1(double x) {
81
82    double w, z, p, q, xn;
83
84    const double DR1 = 5.78318596294678452118E0;
85    const double DR2 = 3.04712623436620863991E1;
86    const double Z1 = 1.46819706421238932572E1;
87    const double Z2 = 4.92184563216946036703E1;
88    const double THPIO4 =  2.35619449019234492885;
89    const double SQ2OPI = 0.79788456080286535588;
90
91    double RP[8] = {
92    -8.99971225705559398224E8,
93    4.52228297998194034323E11,
94    -7.27494245221818276015E13,
95    3.68295732863852883286E15,
96    0.0,
97    0.0,
98    0.0,
99    0.0
100    };
101
102    double RQ[8] = {
103    /* 1.00000000000000000000E0,*/
104    6.20836478118054335476E2,
105    2.56987256757748830383E5,
106    8.35146791431949253037E7,
107    2.21511595479792499675E10,
108    4.74914122079991414898E12,
109    7.84369607876235854894E14,
110    8.95222336184627338078E16,
111    5.32278620332680085395E18,
112    };
113
114    double PP[8] = {
115    7.62125616208173112003E-4,
116    7.31397056940917570436E-2,
117    1.12719608129684925192E0,
118    5.11207951146807644818E0,
119    8.42404590141772420927E0,
120    5.21451598682361504063E0,
121    1.00000000000000000254E0,
122    0.0,
123    };
124    double PQ[8] = {
125    5.71323128072548699714E-4,
126    6.88455908754495404082E-2,
127    1.10514232634061696926E0,
128    5.07386386128601488557E0,
129    8.39985554327604159757E0,
130    5.20982848682361821619E0,
131    9.99999999999999997461E-1,
132    0.0,
133    };
134
135    double QP[8] = {
136    5.10862594750176621635E-2,
137    4.98213872951233449420E0,
138    7.58238284132545283818E1,
139    3.66779609360150777800E2,
140    7.10856304998926107277E2,
141    5.97489612400613639965E2,
142    2.11688757100572135698E2,
143    2.52070205858023719784E1,
144    };
145
146    double QQ[8] = {
147    /* 1.00000000000000000000E0,*/
148    7.42373277035675149943E1,
149    1.05644886038262816351E3,
150    4.98641058337653607651E3,
151    9.56231892404756170795E3,
152    7.99704160447350683650E3,
153    2.82619278517639096600E3,
154    3.36093607810698293419E2,
155    0.0,
156    };
157
158    //const double Z1 = 1.46819706421238932572E1;
159    //const double Z2 = 4.92184563216946036703E1;
160
161    w = x;
162    if( x < 0 )
163            w = -x;
164
165    if( w <= 5.0 )
166        {
167            z = x * x;
168            w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
169            w = w * x * (z - Z1) * (z - Z2);
170            return( w );
171        }
172
173    w = 5.0/x;
174    z = w * w;
175    p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
176    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
177    xn = x - THPIO4;
178
179    double sn, cn;
180    SINCOS(xn, sn, cn);
181    p = p * cn - w * q * sn;
182
183    return( p * SQ2OPI / sqrt(x) );
184}
Note: See TracBrowser for help on using the repository browser.