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


Page generated by Doxygen 1.8.14 for MRPT 1.5.7 Git: 5902e14cc Wed Apr 24 15:04:01 2019 +0200 at lun oct 28 01:39:17 CET 2019