Main MRPT website > C++ reference for MRPT 1.9.9
fresnel.cpp
Go to the documentation of this file.
1 /* +------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | http://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2018, Individual contributors, see AUTHORS file |
6  | See: http://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See details in http://www.mrpt.org/License |
8  +------------------------------------------------------------------------+ */
9 
10 #include "math-precomp.h" // Precompiled headers
11 
12 #include <mrpt/math/fresnel.h>
13 #include <cmath>
14 #include <cfloat>
15 
16 // Based on code (public domain) in http://www.mymathlib.com/functions/
17 
18 static long double Power_Series_S(long double x);
19 static long double xFresnel_Auxiliary_Cosine_Integral(long double x);
20 static long double xFresnel_Auxiliary_Sine_Integral(long double x);
21 static long double Power_Series_C(long double x);
22 static long double xChebyshev_Tn_Series(
23  long double x, const long double a[], int degree);
24 
25 // extern long double xFresnel_Sine_Integral(long double x);
26 // extern long double xFresnel_Cosine_Integral(long double x);
27 
28 // the integral from 0 to x of sqrt(2 / pi) sin(t ^ 2) dt.
29 long double lfresnel_sin_alt(long double x)
30 {
31  long double f, g, x2, s;
32 
33  if (std::abs(x) < 0.5L) return Power_Series_S(x);
34 
35  f = xFresnel_Auxiliary_Cosine_Integral(fabsl(x));
36  g = xFresnel_Auxiliary_Sine_Integral(fabsl(x));
37  x2 = x * x;
38  s = 0.5L - cosl(x2) * f - sinl(x2) * g;
39  return (x < 0.0L) ? -s : s;
40 }
41 
42 // the integral from 0 to x of sqrt(2 / pi) cos(t ^ 2) dt.
43 long double lfresnel_cos_alt(long double x)
44 {
45  long double f, g, x2, c;
46 
47  if (fabsl(x) < 0.5L) return Power_Series_C(x);
48 
49  f = xFresnel_Auxiliary_Cosine_Integral(fabsl(x));
50  g = xFresnel_Auxiliary_Sine_Integral(fabsl(x));
51  x2 = x * x;
52  c = 0.5L + sinl(x2) * f - cosl(x2) * g;
53  return (x < 0.0L) ? -c : c;
54 }
55 
56 long double mrpt::math::lfresnel_sin_integral(long double x) noexcept
57 {
58  long double sqrt_2_o_pi = 7.978845608028653558798921198687637369517e-1L;
59  return lfresnel_sin_alt(x / sqrt_2_o_pi);
60 }
61 
62 long double mrpt::math::lfresnel_cos_integral(long double x) noexcept
63 {
64  long double sqrt_2_o_pi = 7.978845608028653558798921198687637369517e-1L;
65  return lfresnel_cos_alt(x / sqrt_2_o_pi);
66 }
67 
68 double mrpt::math::fresnel_sin_integral(double x) noexcept
69 {
70  return static_cast<double>(
71  lfresnel_sin_integral(static_cast<long double>(x)));
72 }
73 
74 double mrpt::math::fresnel_cos_integral(double x) noexcept
75 {
76  return static_cast<double>(
77  lfresnel_cos_integral(static_cast<long double>(x)));
78 }
79 
80 // Internally Defined Routines //
81 double Fresnel_Sine_Integral(double x);
82 long double xFresnel_Sine_Integral(long double x);
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 // double Fresnel_Sine_Integral( double x ) //
86 // //
87 // Description: //
88 // The Fresnel sine integral, S(x), is the integral with integrand //
89 // sqrt(2/pi) sin(t^2) dt //
90 // where the integral extends from 0 to x. //
91 // //
92 // Arguments: //
93 // double x The argument of the Fresnel sine integral S(). //
94 // //
95 // Return Value: //
96 // The value of the Fresnel sine integral S evaluated at x. //
97 // //
98 // Example: //
99 // double y, x; //
100 // //
101 // ( code to initialize x ) //
102 // //
103 // y = Fresnel_Sine_Integral( x ); //
104 ////////////////////////////////////////////////////////////////////////////////
105 double Fresnel_Sine_Integral(double x)
106 {
107  return (double)xFresnel_Sine_Integral((long double)x);
108 }
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 // long double xFresnel_Sine_Integral( long double x ) //
112 // //
113 // Description: //
114 // The Fresnel sine integral, S(x), is the integral with integrand //
115 // sqrt(2/pi) sin(t^2) dt //
116 // where the integral extends from 0 to x. //
117 // //
118 // Arguments: //
119 // long double x The argument of the Fresnel sine integral S(). //
120 // //
121 // Return Value: //
122 // The value of the Fresnel sine integral S evaluated at x. //
123 // //
124 // Example: //
125 // long double y, x; //
126 // //
127 // ( code to initialize x ) //
128 // //
129 // y = xFresnel_Sine_Integral( x ); //
130 ////////////////////////////////////////////////////////////////////////////////
131 
132 long double xFresnel_Sine_Integral(long double x)
133 {
134  long double f;
135  long double g;
136  long double x2;
137  long double s;
138 
139  if (fabsl(x) < 0.5L) return Power_Series_S(x);
140 
141  f = xFresnel_Auxiliary_Cosine_Integral(fabsl(x));
142  g = xFresnel_Auxiliary_Sine_Integral(fabsl(x));
143  x2 = x * x;
144  s = 0.5L - cosl(x2) * f - sinl(x2) * g;
145  return (x < 0.0L) ? -s : s;
146 }
147 
148 ////////////////////////////////////////////////////////////////////////////////
149 // static long double Power_Series_S( long double x ) //
150 // //
151 // Description: //
152 // The power series representation for the Fresnel sine integral, S(x), //
153 // is //
154 // x^3 sqrt(2/pi) Sum (-x^4)^j / [(4j+3) (2j+1)!] //
155 // where the sum extends over j = 0, ,,,. //
156 // //
157 // Arguments: //
158 // long double x //
159 // The argument of the Fresnel sine integral S(). //
160 // //
161 // Return Value: //
162 // The value of the Fresnel sine integral S evaluated at x. //
163 // //
164 // Example: //
165 // long double y, x; //
166 // //
167 // ( code to initialize x ) //
168 // //
169 // y = Power_Series_S( x ); //
170 ////////////////////////////////////////////////////////////////////////////////
171 
172 static long double Power_Series_S(long double x)
173 {
174  long double x2 = x * x;
175  long double x3 = x * x2;
176  long double x4 = -x2 * x2;
177  long double xn = 1.0L;
178  long double Sn = 1.0L;
179  long double Sm1 = 0.0L;
180  long double term;
181  long double factorial = 1.0L;
182  long double sqrt_2_o_pi = 7.978845608028653558798921198687637369517e-1L;
183  int y = 0;
184 
185  if (x == 0.0L) return 0.0L;
186  Sn /= 3.0L;
187  while (fabsl(Sn - Sm1) > LDBL_EPSILON * fabsl(Sm1))
188  {
189  Sm1 = Sn;
190  y += 1;
191  factorial *= (long double)(y + y);
192  factorial *= (long double)(y + y + 1);
193  xn *= x4;
194  term = xn / factorial;
195  term /= (long double)(y + y + y + y + 3);
196  Sn += term;
197  }
198  return x3 * sqrt_2_o_pi * Sn;
199 }
200 
201 // Internally Defined Routines //
202 double Fresnel_Auxiliary_Sine_Integral(double x);
203 long double xFresnel_Auxiliary_Sine_Integral(long double x);
204 
205 // Internally Defined Constants //
206 static long double const sqrt_2pi = 2.506628274631000502415765284811045253006L;
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 // double Fresnel_Auxiliary_Sine_Integral( double x ) //
210 // //
211 // Description: //
212 // The Fresnel auxiliary sine integral, g(x), is the integral from 0 to //
213 // infinity of the integrand //
214 // sqrt(2/pi) exp(-2xt) sin(t^2) dt //
215 // where x >= 0. //
216 // //
217 // Arguments: //
218 // double x The argument of the Fresnel auxiliary sine integral g() //
219 // where x >= 0. //
220 // //
221 // Return Value: //
222 // The value of the Fresnel auxiliary sine integral g evaluated at //
223 // x >= 0. //
224 // //
225 // Example: //
226 // double y, x; //
227 // //
228 // ( code to initialize x ) //
229 // //
230 // y = Fresnel_Auxiliary_Sine_Integral( x ); //
231 ////////////////////////////////////////////////////////////////////////////////
232 double Fresnel_Auxiliary_Sine_Integral(double x)
233 {
234  return (double)xFresnel_Auxiliary_Sine_Integral((long double)x);
235 }
236 
237 ////////////////////////////////////////////////////////////////////////////////
238 // long double xFresnel_Auxiliary_Sine_Integral( double x ) //
239 // //
240 // Description: //
241 // The Fresnel auxiliary sine integral, g(x), is the integral from 0 to //
242 // infinity of the integrand //
243 // sqrt(2/pi) exp(-2xt) sin(t^2) dt //
244 // where x >= 0. //
245 // //
246 // Arguments: //
247 // long double x The argument of the Fresnel auxiliary sine integral //
248 // g() where x >= 0. //
249 // //
250 // Return Value: //
251 // The value of the Fresnel auxiliary sine integral g evaluated at //
252 // x >= 0. //
253 // //
254 // Example: //
255 // double y, x; //
256 // //
257 // ( code to initialize x ) //
258 // //
259 // y = xFresnel_Auxiliary_Sine_Integral( x ); //
260 ////////////////////////////////////////////////////////////////////////////////
261 
262 static long double sin_Chebyshev_Expansion_0_1(long double x);
263 static long double sin_Chebyshev_Expansion_1_3(long double x);
264 static long double sin_Chebyshev_Expansion_3_5(long double x);
265 static long double sin_Chebyshev_Expansion_5_7(long double x);
266 static long double sin_Asymptotic_Series(long double x);
267 
268 long double xFresnel_Auxiliary_Sine_Integral(long double x)
269 {
270  if (x == 0.0L) return 0.5L;
271  if (x <= 1.0L) return sin_Chebyshev_Expansion_0_1(x);
272  if (x <= 3.0L) return sin_Chebyshev_Expansion_1_3(x);
273  if (x <= 5.0L) return sin_Chebyshev_Expansion_3_5(x);
274  if (x <= 7.0L) return sin_Chebyshev_Expansion_5_7(x);
275  return sin_Asymptotic_Series(x);
276 }
277 
278 ////////////////////////////////////////////////////////////////////////////////
279 // static long double Chebyshev_Expansion_0_1( long double x ) //
280 // //
281 // Description: //
282 // Evaluate the Fresnel auxiliary sine integral, g(x), on the interval //
283 // 0 < x <= 1 using the Chebyshev interpolation formula. //
284 // //
285 // Arguments: //
286 // long double x The argument of the Fresnel auxiliary sine integral //
287 // where 0 < x <= 1. //
288 // //
289 // Return Value: //
290 // The value of the Fresnel auxiliary sine integral g evaluated at //
291 // x where 0 < x <= 1. //
292 // //
293 // Example: //
294 // long double y, x; //
295 // //
296 // ( code to initialize x ) //
297 // //
298 // y = Chebyshev_Expansion_0_1(x); //
299 ////////////////////////////////////////////////////////////////////////////////
300 
301 static long double sin_Chebyshev_Expansion_0_1(long double x)
302 {
303  static long double const c[] = {
304  +2.560134650043040830997e-1L, -1.993005146464943284549e-1L,
305  +4.025503636721387266117e-2L, -4.459600454502960250729e-3L,
306  +6.447097305145147224459e-5L, +7.544218493763717599380e-5L,
307  -1.580422720690700333493e-5L, +1.755845848573471891519e-6L,
308  -9.289769688468301734718e-8L, -5.624033192624251079833e-9L,
309  +1.854740406702369495830e-9L, -2.174644768724492443378e-10L,
310  +1.392899828133395918767e-11L, -6.989216003725983789869e-14L,
311  -9.959396121060010838331e-14L, +1.312085140393647257714e-14L,
312  -9.240470383522792593305e-16L, +2.472168944148817385152e-17L,
313  +2.834615576069400293894e-18L, -4.650983461314449088349e-19L,
314  +3.544083040732391556797e-20L};
315 
316  static const int degree = sizeof(c) / sizeof(long double) - 1;
317  static const long double midpoint = 0.5L;
318  static const long double scale = 0.5L;
319 
320  return xChebyshev_Tn_Series((x - midpoint) / scale, c, degree);
321 }
322 
323 ////////////////////////////////////////////////////////////////////////////////
324 // static long double Chebyshev_Expansion_1_3( long double x ) //
325 // //
326 // Description: //
327 // Evaluate the Fresnel auxiliary sine integral, g(x), on the interval //
328 // 1 < x <= 3 using the Chebyshev interpolation formula. //
329 // //
330 // Arguments: //
331 // long double x The argument of the Fresnel auxiliary sine integral //
332 // where 1 < x <= 3. //
333 // //
334 // Return Value: //
335 // The value of the Fresnel auxiliary sine integral g evaluated at //
336 // x where 1 < x <= 3. //
337 // //
338 // Example: //
339 // long double y, x; //
340 // //
341 // ( code to initialize x ) //
342 // //
343 // y = Chebyshev_Expansion_1_3(x); //
344 ////////////////////////////////////////////////////////////////////////////////
345 
346 static long double sin_Chebyshev_Expansion_1_3(long double x)
347 {
348  static long double const c[] = {
349  +3.470341566046115476477e-2L, -3.855580521778624043304e-2L,
350  +1.420604309383996764083e-2L, -4.037349972538938202143e-3L,
351  +9.292478174580997778194e-4L, -1.742730601244797978044e-4L,
352  +2.563352976720387343201e-5L, -2.498437524746606551732e-6L,
353  -1.334367201897140224779e-8L, +7.436854728157752667212e-8L,
354  -2.059620371321272169176e-8L, +3.753674773239250330547e-9L,
355  -5.052913010605479996432e-10L, +4.580877371233042345794e-11L,
356  -7.664740716178066564952e-13L, -7.200170736686941995387e-13L,
357  +1.812701686438975518372e-13L, -2.799876487275995466163e-14L,
358  +3.048940815174731772007e-15L, -1.936754063718089166725e-16L,
359  -7.653673328908379651914e-18L, +4.534308864750374603371e-18L,
360  -8.011054486030591219007e-19L, +9.374587915222218230337e-20L,
361  -7.144943099280650363024e-21L, +1.105276695821552769144e-22L,
362  +6.989334213887669628647e-23L};
363 
364  static const int degree = sizeof(c) / sizeof(long double) - 1;
365  static const long double midpoint = 2.0L;
366 
367  return xChebyshev_Tn_Series((x - midpoint), c, degree);
368 }
369 
370 ////////////////////////////////////////////////////////////////////////////////
371 // static long double Chebyshev_Expansion_3_5( long double x ) //
372 // //
373 // Description: //
374 // Evaluate the Fresnel auxiliary sine integral, g(x), on the interval //
375 // 3 < x <= 5 using the Chebyshev interpolation formula. //
376 // //
377 // Arguments: //
378 // long double x The argument of the Fresnel auxiliary sine integral //
379 // where 3 < x <= 5. //
380 // //
381 // Return Value: //
382 // The value of the Fresnel auxiliary sine integral g evaluated at //
383 // x where 3 < x <= 5. //
384 // //
385 // Example: //
386 // long double y, x; //
387 // //
388 // ( code to initialize x ) //
389 // //
390 // y = Chebyshev_Expansion_3_5(x); //
391 ////////////////////////////////////////////////////////////////////////////////
392 
393 static long double sin_Chebyshev_Expansion_3_5(long double x)
394 {
395  static long double const c[] = {
396  +3.684922395955255848372e-3L, -2.624595437764014386717e-3L,
397  +6.329162500611499391493e-4L, -1.258275676151483358569e-4L,
398  +2.207375763252044217165e-5L, -3.521929664607266176132e-6L,
399  +5.186211398012883705616e-7L, -7.095056569102400546407e-8L,
400  +9.030550018646936241849e-9L, -1.066057806832232908641e-9L,
401  +1.157128073917012957550e-10L, -1.133877461819345992066e-11L,
402  +9.633572308791154852278e-13L, -6.336675771012312827721e-14L,
403  +1.634407356931822107368e-15L, +3.944542177576016972249e-16L,
404  -9.577486627424256130607e-17L, +1.428772744117447206807e-17L,
405  -1.715342656474756703926e-18L, +1.753564314320837957805e-19L,
406  -1.526125102356904908532e-20L, +1.070275366865736879194e-21L,
407  -4.783978662888842165071e-23L};
408 
409  static const int degree = sizeof(c) / sizeof(long double) - 1;
410  static const long double midpoint = 4.0L;
411 
412  return xChebyshev_Tn_Series((x - midpoint), c, degree);
413 }
414 
415 ////////////////////////////////////////////////////////////////////////////////
416 // static long double Chebyshev_Expansion_5_7( long double x ) //
417 // //
418 // Description: //
419 // Evaluate the Fresnel auxiliary sine integral, g(x), on the interval //
420 // 5 < x <= 7 using the Chebyshev interpolation formula. //
421 // //
422 // Arguments: //
423 // long double x The argument of the Fresnel auxiliary sine integral //
424 // where 5 < x <= 7. //
425 // //
426 // Return Value: //
427 // The value of the Fresnel auxiliary sine integral g evaluated at //
428 // x where 5 < x <= 7. //
429 // //
430 // Example: //
431 // long double y, x; //
432 // //
433 // ( code to initialize x ) //
434 // //
435 // y = Chebyshev_Expansion_5_7(x); //
436 ////////////////////////////////////////////////////////////////////////////////
437 
438 static long double sin_Chebyshev_Expansion_5_7(long double x)
439 {
440  static long double const c[] = {
441  +1.000801217561417083840e-3L, -4.915205279689293180607e-4L,
442  +8.133163567827942356534e-5L, -1.120758739236976144656e-5L,
443  +1.384441872281356422699e-6L, -1.586485067224130537823e-7L,
444  +1.717840749804993618997e-8L, -1.776373217323590289701e-9L,
445  +1.765399783094380160549e-10L, -1.692470022450343343158e-11L,
446  +1.568238301528778401489e-12L, -1.405356860742769958771e-13L,
447  +1.217377701691787512346e-14L, -1.017697418261094517680e-15L,
448  +8.186068056719295045596e-17L, -6.305153620995673221364e-18L,
449  +4.614110100197028845266e-19L, -3.165914620159266813849e-20L,
450  +1.986716456911232767045e-21L, -1.078418278174434671506e-22L,
451  +4.255983404468350776788e-24L};
452 
453  static const int degree = sizeof(c) / sizeof(long double) - 1;
454  static const long double midpoint = 6.0L;
455 
456  return xChebyshev_Tn_Series((x - midpoint), c, degree);
457 }
458 
459 ////////////////////////////////////////////////////////////////////////////////
460 // static long double Asymptotic_Series( long double x ) //
461 // //
462 // Description: //
463 // For a large argument x, the auxiliary Fresnel sine integral, g(x), //
464 // can be expressed as the asymptotic series //
465 // g(x) ~ 1/(x^3 * sqrt(8pi))[1 - 15/4x^4 + 945/16x^8 + ... + //
466 // (4j+1)!!/(-4x^4)^j + ... ] //
467 // //
468 // Arguments: //
469 // long double x The argument of the Fresnel auxiliary sine integral //
470 // where x > 7. //
471 // //
472 // Return Value: //
473 // The value of the Fresnel auxiliary sine integral g evaluated at //
474 // x where x > 7. //
475 // //
476 // Example: //
477 // long double y, x; //
478 // //
479 // ( code to initialize x ) //
480 // //
481 // y = Asymptotic_Series( x ); //
482 ////////////////////////////////////////////////////////////////////////////////
483 #define NUM_ASYMPTOTIC_TERMS 35
484 static long double sin_Asymptotic_Series(long double x)
485 {
486  long double x2 = x * x;
487  long double x4 = -4.0L * x2 * x2;
488  long double xn = 1.0L;
489  long double factorial = 1.0L;
490  long double g = 0.0L;
491  long double term[NUM_ASYMPTOTIC_TERMS + 1];
492  long double epsilon = LDBL_EPSILON / 4.0L;
493  int j = 5;
494  int i = 0;
495 
496  term[0] = 1.0L;
497  term[NUM_ASYMPTOTIC_TERMS] = 0.0L;
498  for (i = 1; i < NUM_ASYMPTOTIC_TERMS; i++)
499  {
500  factorial *= ((long double)j * (long double)(j - 2));
501  xn *= x4;
502  term[i] = factorial / xn;
503  j += 4;
504  if (fabsl(term[i]) >= fabsl(term[i - 1]))
505  {
506  i--;
507  break;
508  }
509  if (fabsl(term[i]) <= epsilon) break;
510  }
511  for (; i >= 0; i--) g += term[i];
512 
513  g /= (x * sqrt_2pi);
514  return g / (x2 + x2);
515 }
516 
517 // Internally Defined Routines //
518 double Fresnel_Auxiliary_Cosine_Integral(double x);
519 long double xFresnel_Auxiliary_Cosine_Integral(long double x);
520 
521 ////////////////////////////////////////////////////////////////////////////////
522 // double Fresnel_Auxiliary_Cosine_Integral( double x ) //
523 // //
524 // Description: //
525 // The Fresnel auxiliary cosine integral, f(x), is the integral from 0 to //
526 // infinity of the integrand //
527 // sqrt(2/pi) exp(-2xt) cos(t^2) dt //
528 // where x >= 0. //
529 // //
530 // Arguments: //
531 // double x The argument of the Fresnel auxiliary cosine integral f() //
532 // where x >= 0. //
533 // //
534 // Return Value: //
535 // The value of the Fresnel auxiliary cosine integral f evaluated at //
536 // x >= 0. //
537 // //
538 // Example: //
539 // double y, x; //
540 // //
541 // ( code to initialize x ) //
542 // //
543 // y = Fresnel_Auxiliary_Cosine_Integral( x ); //
544 ////////////////////////////////////////////////////////////////////////////////
545 double Fresnel_Auxiliary_Cosine_Integral(double x)
546 {
547  return (double)xFresnel_Auxiliary_Cosine_Integral((long double)x);
548 }
549 
550 static long double cos_Chebyshev_Expansion_0_1(long double x);
551 static long double cos_Chebyshev_Expansion_1_3(long double x);
552 static long double cos_Chebyshev_Expansion_3_5(long double x);
553 static long double cos_Chebyshev_Expansion_5_7(long double x);
554 static long double cos_Asymptotic_Series(long double x);
555 
556 ////////////////////////////////////////////////////////////////////////////////
557 // long double xFresnel_Auxiliary_Cosine_Integral( double x ) //
558 // //
559 // Description: //
560 // The Fresnel auxiliary cosine integral, f(x), is the integral from 0 to //
561 // infinity of the integrand //
562 // sqrt(2/pi) exp(-2xt) cos(t^2) dt //
563 // where x >= 0. //
564 // //
565 // Arguments: //
566 // long double x The argument of the Fresnel auxiliary cosine integral //
567 // f() where x >= 0. //
568 // //
569 // Return Value: //
570 // The value of the Fresnel auxiliary cosine integral f evaluated at //
571 // x >= 0. //
572 // //
573 // Example: //
574 // double y, x; //
575 // //
576 // ( code to initialize x ) //
577 // //
578 // y = xFresnel_Auxiliary_Cosine_Integral( x ); //
579 ////////////////////////////////////////////////////////////////////////////////
580 long double xFresnel_Auxiliary_Cosine_Integral(long double x)
581 {
582  if (x == 0.0L) return 0.5L;
583  if (x <= 1.0L) return cos_Chebyshev_Expansion_0_1(x);
584  if (x <= 3.0L) return cos_Chebyshev_Expansion_1_3(x);
585  if (x <= 5.0L) return cos_Chebyshev_Expansion_3_5(x);
586  if (x <= 7.0L) return cos_Chebyshev_Expansion_5_7(x);
587  return cos_Asymptotic_Series(x);
588 }
589 
590 // Externally Defined Routines //
591 extern long double xFresnel_Auxiliary_Cosine_Integral(long double x);
592 extern long double xFresnel_Auxiliary_Sine_Integral(long double x);
593 
594 // Internally Defined Routines //
595 double Fresnel_Cosine_Integral(double x);
596 long double xFresnel_Cosine_Integral(long double x);
597 
598 ////////////////////////////////////////////////////////////////////////////////
599 // double Fresnel_Cosine_Integral( double x ) //
600 // //
601 // Description: //
602 // The Fresnel cosine integral, C(x), is the integral with integrand //
603 // sqrt(2/pi) cos(t^2) dt //
604 // where the integral extends from 0 to x. //
605 // //
606 // Arguments: //
607 // double x The argument of the Fresnel cosine integral C(). //
608 // //
609 // Return Value: //
610 // The value of the Fresnel cosine integral C evaluated at x. //
611 // //
612 // Example: //
613 // double y, x; //
614 // //
615 // ( code to initialize x ) //
616 // //
617 // y = Fresnel_Cosine_Integral( x ); //
618 ////////////////////////////////////////////////////////////////////////////////
619 double Fresnel_Cosine_Integral(double x)
620 {
621  return (double)xFresnel_Cosine_Integral((long double)x);
622 }
623 
624 ////////////////////////////////////////////////////////////////////////////////
625 // long double xFresnel_Cosine_Integral( long double x ) //
626 // //
627 // Description: //
628 // The Fresnel cosine integral, C(x), is the integral with integrand //
629 // sqrt(2/pi) cos(t^2) dt //
630 // where the integral extends from 0 to x. //
631 // //
632 // Arguments: //
633 // long double x The argument of the Fresnel cosine integral C(). //
634 // //
635 // Return Value: //
636 // The value of the Fresnel cosine integral C evaluated at x. //
637 // //
638 // Example: //
639 // long double y, x; //
640 // //
641 // ( code to initialize x ) //
642 // //
643 // y = xFresnel_Cosine_Integral( x ); //
644 ////////////////////////////////////////////////////////////////////////////////
645 
646 long double xFresnel_Cosine_Integral(long double x)
647 {
648  long double f;
649  long double g;
650  long double x2;
651  long double c;
652 
653  if (fabsl(x) < 0.5L) return Power_Series_C(x);
654 
655  f = xFresnel_Auxiliary_Cosine_Integral(fabsl(x));
656  g = xFresnel_Auxiliary_Sine_Integral(fabsl(x));
657  x2 = x * x;
658  c = 0.5L + sinl(x2) * f - cosl(x2) * g;
659  return (x < 0.0L) ? -c : c;
660 }
661 
662 ////////////////////////////////////////////////////////////////////////////////
663 // static long double Power_Series_C( long double x ) //
664 // //
665 // Description: //
666 // The power series representation for the Fresnel cosine integral, C(x), //
667 // is //
668 // x sqrt(2/pi) Sum (-x^4)^j / [(4j+1) (2j)!] //
669 // where the sum extends over j = 0, ,,,. //
670 // //
671 // Arguments: //
672 // long double x //
673 // The argument of the Fresnel cosine integral C(). //
674 // //
675 // Return Value: //
676 // The value of the Fresnel cosine integral C evaluated at x. //
677 // //
678 // Example: //
679 // long double y, x; //
680 // //
681 // ( code to initialize x ) //
682 // //
683 // y = Power_Series_C( x ); //
684 ////////////////////////////////////////////////////////////////////////////////
685 
686 static long double Power_Series_C(long double x)
687 {
688  long double x2 = x * x;
689  // long double x3 = x * x2;
690  long double x4 = -x2 * x2;
691  long double xn = 1.0L;
692  long double Sn = 1.0L;
693  long double Sm1 = 0.0L;
694  long double term;
695  long double factorial = 1.0L;
696  long double sqrt_2_o_pi = 7.978845608028653558798921198687637369517e-1L;
697  int y = 0;
698 
699  if (x == 0.0L) return 0.0L;
700  while (fabsl(Sn - Sm1) > LDBL_EPSILON * fabsl(Sm1))
701  {
702  Sm1 = Sn;
703  y += 1;
704  factorial *= (long double)(y + y);
705  factorial *= (long double)(y + y - 1);
706  xn *= x4;
707  term = xn / factorial;
708  term /= (long double)(y + y + y + y + 1);
709  Sn += term;
710  }
711  return x * sqrt_2_o_pi * Sn;
712 }
713 
714 ////////////////////////////////////////////////////////////////////////////////
715 // File: xchebyshev_Tn_series.c //
716 // Routine(s): //
717 // xChebyshev_Tn_Series //
718 ////////////////////////////////////////////////////////////////////////////////
719 ////////////////////////////////////////////////////////////////////////////////
720 // long double xChebyshev_Tn_Series(long double x, long double a[],int degree)//
721 // //
722 // Description: //
723 // This routine uses Clenshaw's recursion algorithm to evaluate a given //
724 // polynomial p(x) expressed as a linear combination of Chebyshev //
725 // polynomials of the first kind, Tn, at a point x, //
726 // p(x) = a[0] + a[1]*T[1](x) + a[2]*T[2](x) + ... + a[deg]*T[deg](x). //
727 // //
728 // Clenshaw's recursion formula applied to Chebyshev polynomials of the //
729 // first kind is: //
730 // Set y[degree + 2] = 0, y[degree + 1] = 0, then for k = degree, ..., 1 //
731 // set y[k] = 2 * x * y[k+1] - y[k+2] + a[k]. Finally //
732 // set y[0] = x * y[1] - y[2] + a[0]. Then p(x) = y[0]. //
733 // //
734 // Arguments: //
735 // long double x //
736 // The point at which to evaluate the polynomial. //
737 // long double a[] //
738 // The coefficients of the expansion in terms of Chebyshev polynomials,//
739 // i.e. a[k] is the coefficient of T[k](x). Note that in the calling //
740 // routine a must be defined double a[N] where N >= degree + 1. //
741 // int degree //
742 // The degree of the polynomial p(x). //
743 // //
744 // Return Value: //
745 // The value of the polynomial at x. //
746 // If degree is negative, then 0.0 is returned. //
747 // //
748 // Example: //
749 // long double x, a[N], p; //
750 // int deg = N - 1; //
751 // //
752 // ( code to initialize x, and a[i] i = 0, ... , a[deg] ) //
753 // //
754 // p = xChebyshev_Tn_Series(x, a, deg); //
755 ////////////////////////////////////////////////////////////////////////////////
756 static long double xChebyshev_Tn_Series(
757  long double x, const long double a[], int degree)
758 {
759  long double yp2 = 0.0L;
760  long double yp1 = 0.0L;
761  long double y = 0.0L;
762  long double two_x = x + x;
763  int k;
764 
765  // Check that degree >= 0. If not, then return 0. //
766 
767  if (degree < 0) return 0.0L;
768 
769  // Apply Clenshaw's recursion save the last iteration. //
770 
771  for (k = degree; k >= 1; k--, yp2 = yp1, yp1 = y)
772  y = two_x * yp1 - yp2 + a[k];
773 
774  // Now apply the last iteration and return the result. //
775 
776  return x * yp1 - yp2 + a[0];
777 }
778 
779 ////////////////////////////////////////////////////////////////////////////////
780 // static long double Chebyshev_Expansion_0_1( long double x ) //
781 // //
782 // Description: //
783 // Evaluate the Fresnel auxiliary cosine integral, f(x), on the interval //
784 // 0 < x <= 1 using the Chebyshev interpolation formula. //
785 // //
786 // Arguments: //
787 // long double x The argument of the Fresnel auxiliary cosine integral //
788 // where 0 < x <= 1. //
789 // //
790 // Return Value: //
791 // The value of the Fresnel auxiliary cosine integral f evaluated at //
792 // x where 0 < x <= 1. //
793 // //
794 // Example: //
795 // long double y, x; //
796 // //
797 // ( code to initialize x ) //
798 // //
799 // y = Chebyshev_Expansion_0_1(x); //
800 ////////////////////////////////////////////////////////////////////////////////
801 
802 static long double cos_Chebyshev_Expansion_0_1(long double x)
803 {
804  static long double const c[] = {
805  +4.200987560240514577713e-1L, -9.358785913634965235904e-2L,
806  -7.642539415723373644927e-3L, +4.958117751796130135544e-3L,
807  -9.750236036106120253456e-4L, +1.075201474958704192865e-4L,
808  -4.415344769301324238886e-6L, -7.861633919783064216022e-7L,
809  +1.919240966215861471754e-7L, -2.175775608982741065385e-8L,
810  +1.296559541430849437217e-9L, +2.207205095025162212169e-11L,
811  -1.479219615873704298874e-11L, +1.821350127295808288614e-12L,
812  -1.228919312990171362342e-13L, +2.227139250593818235212e-15L,
813  +5.734729405928016301596e-16L, -8.284965573075354177016e-17L,
814  +6.067422701530157308321e-18L, -1.994908519477689596319e-19L,
815  -1.173365630675305693390e-20L};
816 
817  static const int degree = sizeof(c) / sizeof(long double) - 1;
818  static const long double midpoint = 0.5L;
819  static const long double scale = 0.5L;
820 
821  return xChebyshev_Tn_Series((x - midpoint) / scale, c, degree);
822 }
823 
824 ////////////////////////////////////////////////////////////////////////////////
825 // static long double Chebyshev_Expansion_1_3( long double x ) //
826 // //
827 // Description: //
828 // Evaluate the Fresnel auxiliary cosine integral, f(x), on the interval //
829 // 1 < x <= 3 using the Chebyshev interpolation formula. //
830 // //
831 // Arguments: //
832 // long double x The argument of the Fresnel auxiliary cosine integral //
833 // where 1 < x <= 3. //
834 // //
835 // Return Value: //
836 // The value of the Fresnel auxiliary cosine integral f evaluated at //
837 // x where 1 < x <= 3. //
838 // //
839 // Example: //
840 // long double y, x; //
841 // //
842 // ( code to initialize x ) //
843 // //
844 // y = Chebyshev_Expansion_1_3(x); //
845 ////////////////////////////////////////////////////////////////////////////////
846 
847 static long double cos_Chebyshev_Expansion_1_3(long double x)
848 {
849  static long double const c[] = {
850  +2.098677278318224971989e-1L, -9.314234883154103266195e-2L,
851  +1.739905936938124979297e-2L, -2.454274824644285136137e-3L,
852  +1.589872606981337312438e-4L, +4.203943842506079780413e-5L,
853  -2.018022256093216535093e-5L, +5.125709636776428285284e-6L,
854  -9.601813551752718650057e-7L, +1.373989484857155846826e-7L,
855  -1.348105546577211255591e-8L, +2.745868700337953872632e-10L,
856  +2.401655517097260106976e-10L, -6.678059547527685587692e-11L,
857  +1.140562171732840809159e-11L, -1.401526517205212219089e-12L,
858  +1.105498827380224475667e-13L, +2.040731455126809208066e-16L,
859  -1.946040679213045143184e-15L, +4.151821375667161733612e-16L,
860  -5.642257647205149369594e-17L, +5.266176626521504829010e-18L,
861  -2.299025577897146333791e-19L, -2.952226367506641078731e-20L,
862  +8.760405943193778149078e-21L};
863 
864  static const int degree = sizeof(c) / sizeof(long double) - 1;
865  static const long double midpoint = 2.0L;
866 
867  return xChebyshev_Tn_Series((x - midpoint), c, degree);
868 }
869 
870 ////////////////////////////////////////////////////////////////////////////////
871 // static long double Chebyshev_Expansion_3_5( long double x ) //
872 // //
873 // Description: //
874 // Evaluate the Fresnel auxiliary cosine integral, g(x), on the interval //
875 // 3 < x <= 5 using the Chebyshev interpolation formula. //
876 // //
877 // Arguments: //
878 // long double x The argument of the Fresnel auxiliary cosine integral //
879 // where 3 < x <= 5. //
880 // //
881 // Return Value: //
882 // The value of the Fresnel auxiliary cosine integral f evaluated at //
883 // x where 3 < x <= 5. //
884 // //
885 // Example: //
886 // long double y, x; //
887 // //
888 // ( code to initialize x ) //
889 // //
890 // y = Chebyshev_Expansion_3_5(x); //
891 ////////////////////////////////////////////////////////////////////////////////
892 
893 static long double cos_Chebyshev_Expansion_3_5(long double x)
894 {
895  static long double const c[] = {
896  +1.025703371090289562388e-1L, -2.569833023232301400495e-2L,
897  +3.160592981728234288078e-3L, -3.776110718882714758799e-4L,
898  +4.325593433537248833341e-5L, -4.668447489229591855730e-6L,
899  +4.619254757356785108280e-7L, -3.970436510433553795244e-8L,
900  +2.535664754977344448598e-9L, -2.108170964644819803367e-11L,
901  -2.959172018518707683013e-11L, +6.727219944906606516055e-12L,
902  -1.062829587519902899001e-12L, +1.402071724705287701110e-13L,
903  -1.619154679722651005075e-14L, +1.651319588396970446858e-15L,
904  -1.461704569438083772889e-16L, +1.053521559559583268504e-17L,
905  -4.760946403462515858756e-19L, -1.803784084922403924313e-20L,
906  +7.873130866418738207547e-21L};
907 
908  static const int degree = sizeof(c) / sizeof(long double) - 1;
909  static const long double midpoint = 4.0L;
910 
911  return xChebyshev_Tn_Series((x - midpoint), c, degree);
912 }
913 
914 ////////////////////////////////////////////////////////////////////////////////
915 // static long double Chebyshev_Expansion_5_7( long double x ) //
916 // //
917 // Description: //
918 // Evaluate the Fresnel auxiliary cosine integral, g(x), on the interval //
919 // 5 < x <= 7 using the Chebyshev interpolation formula. //
920 // //
921 // Arguments: //
922 // long double x The argument of the Fresnel auxiliary cosine integral //
923 // where 5 < x <= 7. //
924 // //
925 // Return Value: //
926 // The value of the Fresnel auxiliary cosine integral f evaluated at //
927 // x where 5 < x <= 7. //
928 // //
929 // Example: //
930 // long double y, x; //
931 // //
932 // ( code to initialize x ) //
933 // //
934 // y = Chebyshev_Expansion_5_7(x); //
935 ////////////////////////////////////////////////////////////////////////////////
936 
937 static long double cos_Chebyshev_Expansion_5_7(long double x)
938 {
939  static long double const c[] = {
940  +6.738667333400589274018e-2L, -1.128146832637904868638e-2L,
941  +9.408843234170404670278e-4L, -7.800074103496165011747e-5L,
942  +6.409101169623350885527e-6L, -5.201350558247239981834e-7L,
943  +4.151668914650221476906e-8L, -3.242202015335530552721e-9L,
944  +2.460339340900396789789e-10L, -1.796823324763304661865e-11L,
945  +1.244108496436438952425e-12L, -7.950417122987063540635e-14L,
946  +4.419142625999150971878e-15L, -1.759082736751040110146e-16L,
947  -1.307443936270786700760e-18L, +1.362484141039320395814e-18L,
948  -2.055236564763877250559e-19L, +2.329142055084791308691e-20L,
949  -2.282438671525884861970e-21L};
950 
951  static const int degree = sizeof(c) / sizeof(long double) - 1;
952  static const long double midpoint = 6.0L;
953 
954  return xChebyshev_Tn_Series((x - midpoint), c, degree);
955 }
956 
957 ////////////////////////////////////////////////////////////////////////////////
958 // static long double Asymptotic_Series( long double x ) //
959 // //
960 // Description: //
961 // For a large argument x, the auxiliary Fresnel cosine integral, f(x), //
962 // can be expressed as the asymptotic series //
963 // f(x) ~ 1/(x*sqrt(2pi))[1 - 3/4x^4 + 105/16x^8 + ... + //
964 // (4j-1)!!/(-4x^4)^j + ... ] //
965 // //
966 // Arguments: //
967 // long double x The argument of the Fresnel auxiliary cosine integral //
968 // where x > 7. //
969 // //
970 // Return Value: //
971 // The value of the Fresnel auxiliary cosine integral f evaluated at //
972 // x where x > 7. //
973 // //
974 // Example: //
975 // long double y, x; //
976 // //
977 // ( code to initialize x ) //
978 // //
979 // y = Asymptotic_Series( x ); //
980 ////////////////////////////////////////////////////////////////////////////////
981 #define NUM_ASYMPTOTIC_TERMS 35
982 static long double cos_Asymptotic_Series(long double x)
983 {
984  long double x2 = x * x;
985  long double x4 = -4.0L * x2 * x2;
986  long double xn = 1.0L;
987  long double factorial = 1.0L;
988  long double f = 0.0L;
989  long double term[NUM_ASYMPTOTIC_TERMS + 1];
990  long double epsilon = LDBL_EPSILON / 4.0L;
991  int j = 3;
992  int i = 0;
993 
994  term[0] = 1.0L;
995  term[NUM_ASYMPTOTIC_TERMS] = 0.0L;
996  for (i = 1; i < NUM_ASYMPTOTIC_TERMS; i++)
997  {
998  factorial *= ((long double)j * (long double)(j - 2));
999  xn *= x4;
1000  term[i] = factorial / xn;
1001  j += 4;
1002  if (fabsl(term[i]) >= fabsl(term[i - 1]))
1003  {
1004  i--;
1005  break;
1006  }
1007  if (fabsl(term[i]) <= epsilon) break;
1008  }
1009 
1010  for (; i >= 0; i--) f += term[i];
1011 
1012  return f / (x * sqrt_2pi);
1013 }
cos_Chebyshev_Expansion_1_3
static long double cos_Chebyshev_Expansion_1_3(long double x)
Definition: fresnel.cpp:847
Fresnel_Sine_Integral
double Fresnel_Sine_Integral(double x)
Definition: fresnel.cpp:105
cos_Asymptotic_Series
static long double cos_Asymptotic_Series(long double x)
Definition: fresnel.cpp:982
s
GLdouble s
Definition: glext.h:3676
mrpt::math::factorial
double factorial(unsigned int n)
Computes the factorial of an integer number and returns it as a double value (internally it uses loga...
Definition: math.cpp:68
c
const GLubyte * c
Definition: glext.h:6313
sin_Chebyshev_Expansion_5_7
static long double sin_Chebyshev_Expansion_5_7(long double x)
Definition: fresnel.cpp:438
math-precomp.h
sin_Chebyshev_Expansion_3_5
static long double sin_Chebyshev_Expansion_3_5(long double x)
Definition: fresnel.cpp:393
fresnel.h
sqrt_2pi
static long const double sqrt_2pi
Definition: fresnel.cpp:206
cos_Chebyshev_Expansion_0_1
static long double cos_Chebyshev_Expansion_0_1(long double x)
Definition: fresnel.cpp:802
lfresnel_cos_alt
long double lfresnel_cos_alt(long double x)
Definition: fresnel.cpp:43
Power_Series_S
static long double Power_Series_S(long double x)
Definition: fresnel.cpp:172
g
GLubyte g
Definition: glext.h:6279
mrpt::math::lfresnel_sin_integral
long double lfresnel_sin_integral(long double x) noexcept
long double version of fresnel_sin_integral
Definition: fresnel.cpp:56
cos_Chebyshev_Expansion_5_7
static long double cos_Chebyshev_Expansion_5_7(long double x)
Definition: fresnel.cpp:937
sin_Asymptotic_Series
static long double sin_Asymptotic_Series(long double x)
Definition: fresnel.cpp:484
Fresnel_Cosine_Integral
double Fresnel_Cosine_Integral(double x)
Definition: fresnel.cpp:619
mrpt::math::lfresnel_cos_integral
long double lfresnel_cos_integral(long double x) noexcept
long double version of fresnel_cos_integral
Definition: fresnel.cpp:62
cos_Chebyshev_Expansion_3_5
static long double cos_Chebyshev_Expansion_3_5(long double x)
Definition: fresnel.cpp:893
lfresnel_sin_alt
long double lfresnel_sin_alt(long double x)
Definition: fresnel.cpp:29
mrpt::math::fresnel_sin_integral
double fresnel_sin_integral(double x) noexcept
Evaluates the integral from 0 to x of sqrt(2/pi) sin(t^2) dt.
Definition: fresnel.cpp:68
xFresnel_Sine_Integral
long double xFresnel_Sine_Integral(long double x)
Definition: fresnel.cpp:132
NUM_ASYMPTOTIC_TERMS
#define NUM_ASYMPTOTIC_TERMS
Definition: fresnel.cpp:981
Power_Series_C
static long double Power_Series_C(long double x)
Definition: fresnel.cpp:686
scale
GLenum GLenum GLenum GLenum GLenum scale
Definition: glext.h:6503
xFresnel_Cosine_Integral
long double xFresnel_Cosine_Integral(long double x)
Definition: fresnel.cpp:646
sin_Chebyshev_Expansion_1_3
static long double sin_Chebyshev_Expansion_1_3(long double x)
Definition: fresnel.cpp:346
xFresnel_Auxiliary_Cosine_Integral
static long double xFresnel_Auxiliary_Cosine_Integral(long double x)
Definition: fresnel.cpp:580
Fresnel_Auxiliary_Cosine_Integral
double Fresnel_Auxiliary_Cosine_Integral(double x)
Definition: fresnel.cpp:545
y
GLenum GLint GLint y
Definition: glext.h:3538
mrpt::math::fresnel_cos_integral
double fresnel_cos_integral(double x) noexcept
Evaluates the integral from 0 to x of sqrt(2/pi) cos(t^2) dt.
Definition: fresnel.cpp:74
x
GLenum GLint x
Definition: glext.h:3538
xChebyshev_Tn_Series
static long double xChebyshev_Tn_Series(long double x, const long double a[], int degree)
Definition: fresnel.cpp:756
Fresnel_Auxiliary_Sine_Integral
double Fresnel_Auxiliary_Sine_Integral(double x)
Definition: fresnel.cpp:232
xFresnel_Auxiliary_Sine_Integral
static long double xFresnel_Auxiliary_Sine_Integral(long double x)
Definition: fresnel.cpp:268
a
GLubyte GLubyte GLubyte a
Definition: glext.h:6279
sin_Chebyshev_Expansion_0_1
static long double sin_Chebyshev_Expansion_0_1(long double x)
Definition: fresnel.cpp:301


Page generated by Doxygen 1.8.17 for MRPT 1.9.9 Git: ad3a9d8ae Tue May 1 23:10:22 2018 -0700 at mié 12 jul 2023 10:03:34 CEST