fresnl.c
Go to the documentation of this file.
1 /* fresnl.c
2  *
3  * Fresnel integral
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, S, C;
10  * void fresnl();
11  *
12  * fresnl( x, _&S, _&C );
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Evaluates the Fresnel integrals
18  *
19  * x
20  * -
21  * | |
22  * C(x) = | cos(pi/2 t**2) dt,
23  * | |
24  * -
25  * 0
26  *
27  * x
28  * -
29  * | |
30  * S(x) = | sin(pi/2 t**2) dt.
31  * | |
32  * -
33  * 0
34  *
35  *
36  * The integrals are evaluated by a power series for x < 1.
37  * For x >= 1 auxiliary functions f(x) and g(x) are employed
38  * such that
39  *
40  * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
41  * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
42  *
43  *
44  *
45  * ACCURACY:
46  *
47  * Relative error.
48  *
49  * Arithmetic function domain # trials peak rms
50  * IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16
51  * IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16
52  */
53 
54 /*
55  * Cephes Math Library Release 2.1: January, 1989
56  * Copyright 1984, 1987, 1989 by Stephen L. Moshier
57  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
58  */
59 
60 #include "mconf.h"
61 
62 /* S(x) for small x */
63 static double sn[6] = {
64  -2.99181919401019853726E3,
65  7.08840045257738576863E5,
66  -6.29741486205862506537E7,
67  2.54890880573376359104E9,
68  -4.42979518059697779103E10,
69  3.18016297876567817986E11,
70 };
71 
72 static double sd[6] = {
73  /* 1.00000000000000000000E0, */
74  2.81376268889994315696E2,
75  4.55847810806532581675E4,
76  5.17343888770096400730E6,
77  4.19320245898111231129E8,
78  2.24411795645340920940E10,
79  6.07366389490084639049E11,
80 };
81 
82 /* C(x) for small x */
83 static double cn[6] = {
84  -4.98843114573573548651E-8,
85  9.50428062829859605134E-6,
86  -6.45191435683965050962E-4,
87  1.88843319396703850064E-2,
88  -2.05525900955013891793E-1,
89  9.99999999999999998822E-1,
90 };
91 
92 static double cd[7] = {
93  3.99982968972495980367E-12,
94  9.15439215774657478799E-10,
95  1.25001862479598821474E-7,
96  1.22262789024179030997E-5,
97  8.68029542941784300606E-4,
98  4.12142090722199792936E-2,
99  1.00000000000000000118E0,
100 };
101 
102 /* Auxiliary function f(x) */
103 static double fn[10] = {
104  4.21543555043677546506E-1,
105  1.43407919780758885261E-1,
106  1.15220955073585758835E-2,
107  3.45017939782574027900E-4,
108  4.63613749287867322088E-6,
109  3.05568983790257605827E-8,
110  1.02304514164907233465E-10,
111  1.72010743268161828879E-13,
112  1.34283276233062758925E-16,
113  3.76329711269987889006E-20,
114 };
115 
116 static double fd[10] = {
117  /* 1.00000000000000000000E0, */
118  7.51586398353378947175E-1,
119  1.16888925859191382142E-1,
120  6.44051526508858611005E-3,
121  1.55934409164153020873E-4,
122  1.84627567348930545870E-6,
123  1.12699224763999035261E-8,
124  3.60140029589371370404E-11,
125  5.88754533621578410010E-14,
126  4.52001434074129701496E-17,
127  1.25443237090011264384E-20,
128 };
129 
130 /* Auxiliary function g(x) */
131 static double gn[11] = {
132  5.04442073643383265887E-1,
133  1.97102833525523411709E-1,
134  1.87648584092575249293E-2,
135  6.84079380915393090172E-4,
136  1.15138826111884280931E-5,
137  9.82852443688422223854E-8,
138  4.45344415861750144738E-10,
139  1.08268041139020870318E-12,
140  1.37555460633261799868E-15,
141  8.36354435630677421531E-19,
142  1.86958710162783235106E-22,
143 };
144 
145 static double gd[11] = {
146  /* 1.00000000000000000000E0, */
147  1.47495759925128324529E0,
148  3.37748989120019970451E-1,
149  2.53603741420338795122E-2,
150  8.14679107184306179049E-4,
151  1.27545075667729118702E-5,
152  1.04314589657571990585E-7,
153  4.60680728146520428211E-10,
154  1.10273215066240270757E-12,
155  1.38796531259578871258E-15,
156  8.39158816283118707363E-19,
157  1.86958710162783236342E-22,
158 };
159 
160 extern double MACHEP;
161 
162 int fresnl(double xxa, double *ssa, double *cca)
163 {
164  double f, g, cc, ss, c, s, t, u;
165  double x, x2;
166 
167  if (cephes_isinf(xxa)) {
168  cc = 0.5;
169  ss = 0.5;
170  goto done;
171  }
172 
173  x = fabs(xxa);
174  x2 = x * x;
175  if (x2 < 2.5625) {
176  t = x2 * x2;
177  ss = x * x2 * polevl(t, sn, 5) / p1evl(t, sd, 6);
178  cc = x * polevl(t, cn, 5) / polevl(t, cd, 6);
179  goto done;
180  }
181 
182  if (x > 36974.0) {
183  /*
184  * http://functions.wolfram.com/GammaBetaErf/FresnelC/06/02/
185  * http://functions.wolfram.com/GammaBetaErf/FresnelS/06/02/
186  */
187  cc = 0.5 + 1/(M_PI*x) * sin(M_PI*x*x/2);
188  ss = 0.5 - 1/(M_PI*x) * cos(M_PI*x*x/2);
189  goto done;
190  }
191 
192 
193  /* Asymptotic power series auxiliary functions
194  * for large argument
195  */
196  x2 = x * x;
197  t = M_PI * x2;
198  u = 1.0 / (t * t);
199  t = 1.0 / t;
200  f = 1.0 - u * polevl(u, fn, 9) / p1evl(u, fd, 10);
201  g = t * polevl(u, gn, 10) / p1evl(u, gd, 11);
202 
203  t = M_PI_2 * x2;
204  c = cos(t);
205  s = sin(t);
206  t = M_PI * x;
207  cc = 0.5 + (f * s - g * c) / t;
208  ss = 0.5 - (f * c + g * s) / t;
209 
210  done:
211  if (xxa < 0.0) {
212  cc = -cc;
213  ss = -ss;
214  }
215 
216  *cca = cc;
217  *ssa = ss;
218  return (0);
219 }
gn
static double gn[11]
Definition: fresnl.c:131
MACHEP
double MACHEP
Definition: const.c:54
s
RealScalar s
Definition: level1_cplx_impl.h:126
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
fn
static double fn[10]
Definition: fresnl.c:103
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
fresnl
int fresnl(double xxa, double *ssa, double *cca)
Definition: fresnl.c:162
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
cd
static double cd[7]
Definition: fresnl.c:92
cephes_isinf
#define cephes_isinf(x)
Definition: mconf.h:100
sd
static double sd[6]
Definition: fresnl.c:72
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
ss
static std::stringstream ss
Definition: testBTree.cpp:31
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
cn
static double cn[6]
Definition: fresnl.c:83
g
void g(const string &key, int i)
Definition: testBTree.cpp:41
M_PI_2
#define M_PI_2
Definition: mconf.h:118
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
sn
static double sn[6]
Definition: fresnl.c:63
mconf.h
fd
static double fd[10]
Definition: fresnl.c:116
M_PI
#define M_PI
Definition: mconf.h:117
gd
static double gd[11]
Definition: fresnl.c:145
align_3::t
Point2 t(10, 10)
x2
Pose3 x2(Rot3::Ypr(0.0, 0.0, 0.0), l2)


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