j1.c
Go to the documentation of this file.
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  * IEEE 0, 30 30000 2.6e-16 1.1e-16
32  *
33  *
34  */
35  /* y1.c
36  *
37  * Bessel function of second kind of order one
38  *
39  *
40  *
41  * SYNOPSIS:
42  *
43  * double x, y, y1();
44  *
45  * y = y1( x );
46  *
47  *
48  *
49  * DESCRIPTION:
50  *
51  * Returns Bessel function of the second kind of order one
52  * of the argument.
53  *
54  * The domain is divided into the intervals [0, 8] and
55  * (8, infinity). In the first interval a 25 term Chebyshev
56  * expansion is used, and a call to j1() is required.
57  * In the second, the asymptotic trigonometric representation
58  * is employed using two rational functions of degree 5/5.
59  *
60  *
61  *
62  * ACCURACY:
63  *
64  * Absolute error:
65  * arithmetic domain # trials peak rms
66  * IEEE 0, 30 30000 1.0e-15 1.3e-16
67  *
68  * (error criterion relative when |y1| > 1).
69  *
70  */
71 
72 
73 /*
74  * Cephes Math Library Release 2.8: June, 2000
75  * Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
76  */
77 
78 /*
79  * #define PIO4 .78539816339744830962
80  * #define THPIO4 2.35619449019234492885
81  * #define SQ2OPI .79788456080286535588
82  */
83 
84 #include "mconf.h"
85 
86 static double RP[4] = {
87  -8.99971225705559398224E8,
88  4.52228297998194034323E11,
89  -7.27494245221818276015E13,
90  3.68295732863852883286E15,
91 };
92 
93 static double RQ[8] = {
94  /* 1.00000000000000000000E0, */
95  6.20836478118054335476E2,
96  2.56987256757748830383E5,
97  8.35146791431949253037E7,
98  2.21511595479792499675E10,
99  4.74914122079991414898E12,
100  7.84369607876235854894E14,
101  8.95222336184627338078E16,
102  5.32278620332680085395E18,
103 };
104 
105 static double PP[7] = {
106  7.62125616208173112003E-4,
107  7.31397056940917570436E-2,
108  1.12719608129684925192E0,
109  5.11207951146807644818E0,
110  8.42404590141772420927E0,
111  5.21451598682361504063E0,
112  1.00000000000000000254E0,
113 };
114 
115 static double PQ[7] = {
116  5.71323128072548699714E-4,
117  6.88455908754495404082E-2,
118  1.10514232634061696926E0,
119  5.07386386128601488557E0,
120  8.39985554327604159757E0,
121  5.20982848682361821619E0,
122  9.99999999999999997461E-1,
123 };
124 
125 static double QP[8] = {
126  5.10862594750176621635E-2,
127  4.98213872951233449420E0,
128  7.58238284132545283818E1,
129  3.66779609360150777800E2,
130  7.10856304998926107277E2,
131  5.97489612400613639965E2,
132  2.11688757100572135698E2,
133  2.52070205858023719784E1,
134 };
135 
136 static double QQ[7] = {
137  /* 1.00000000000000000000E0, */
138  7.42373277035675149943E1,
139  1.05644886038262816351E3,
140  4.98641058337653607651E3,
141  9.56231892404756170795E3,
142  7.99704160447350683650E3,
143  2.82619278517639096600E3,
144  3.36093607810698293419E2,
145 };
146 
147 static double YP[6] = {
148  1.26320474790178026440E9,
149  -6.47355876379160291031E11,
150  1.14509511541823727583E14,
151  -8.12770255501325109621E15,
152  2.02439475713594898196E17,
153  -7.78877196265950026825E17,
154 };
155 
156 static double YQ[8] = {
157  /* 1.00000000000000000000E0, */
158  5.94301592346128195359E2,
159  2.35564092943068577943E5,
160  7.34811944459721705660E7,
161  1.87601316108706159478E10,
162  3.88231277496238566008E12,
163  6.20557727146953693363E14,
164  6.87141087355300489866E16,
165  3.97270608116560655612E18,
166 };
167 
168 
169 static double Z1 = 1.46819706421238932572E1;
170 static double Z2 = 4.92184563216946036703E1;
171 
172 extern double THPIO4, SQ2OPI;
173 
174 double j1(double x)
175 {
176  double w, z, p, q, xn;
177 
178  w = x;
179  if (x < 0)
180  return -j1(-x);
181 
182  if (w <= 5.0) {
183  z = x * x;
184  w = polevl(z, RP, 3) / p1evl(z, RQ, 8);
185  w = w * x * (z - Z1) * (z - Z2);
186  return (w);
187  }
188 
189  w = 5.0 / x;
190  z = w * w;
191  p = polevl(z, PP, 6) / polevl(z, PQ, 6);
192  q = polevl(z, QP, 7) / p1evl(z, QQ, 7);
193  xn = x - THPIO4;
194  p = p * cos(xn) - w * q * sin(xn);
195  return (p * SQ2OPI / sqrt(x));
196 }
197 
198 
199 double y1(double x)
200 {
201  double w, z, p, q, xn;
202 
203  if (x <= 5.0) {
204  if (x == 0.0) {
206  return -INFINITY;
207  }
208  else if (x <= 0.0) {
209  sf_error("y1", SF_ERROR_DOMAIN, NULL);
210  return NAN;
211  }
212  z = x * x;
213  w = x * (polevl(z, YP, 5) / p1evl(z, YQ, 8));
214  w += M_2_PI * (j1(x) * log(x) - 1.0 / x);
215  return (w);
216  }
217 
218  w = 5.0 / x;
219  z = w * w;
220  p = polevl(z, PP, 6) / polevl(z, PQ, 6);
221  q = polevl(z, QP, 7) / p1evl(z, QQ, 7);
222  xn = x - THPIO4;
223  p = p * sin(xn) + w * q * cos(xn);
224  return (p * SQ2OPI / sqrt(x));
225 }
QP
static double QP[8]
Definition: j1.c:125
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
Z1
static double Z1
Definition: j1.c:169
THPIO4
double THPIO4
Definition: const.c:67
PQ
static double PQ[7]
Definition: j1.c:115
RQ
static double RQ[8]
Definition: j1.c:93
M_2_PI
#define M_2_PI
Definition: mconf.h:121
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Z2
Cyclic< 2 > Z2
Definition: testCyclic.cpp:26
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
RP
static double RP[4]
Definition: j1.c:86
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
SQ2OPI
double SQ2OPI
Definition: j1.c:172
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
y1
double y1(double x)
Definition: j1.c:199
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
YQ
static double YQ[8]
Definition: j1.c:156
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
PP
static double PP[7]
Definition: j1.c:105
YP
static double YP[6]
Definition: j1.c:147
mconf.h
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
j1
double j1(double x)
Definition: j1.c:174
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
QQ
static double QQ[7]
Definition: j1.c:136
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:55