Fresnel integrals

Calculates FresnelC, FresnelS integrals and their variations (C_1, C_2 and S_1 and S_2)

Description

FRESNEL(X) calculates the values of the Fresnel integrals for real values of vector X,
   i.e.
       C = \int_0^x cos(pi*t^2/2) dt, (0a)
       S = \int_0^x sin(pi*t^2/2) dt (0b)
Also, it evaluates the following variations of the Fresnel integrals
       C1 = \sqrt(2/pi) \int_0^x cos(t^2) dt, (1a)
       S1 = \sqrt(2/pi) \int_0^x sin(t^2) dt (1b)
   and
       C2 = \sqrt(1/2/pi) \int_0^x cos(t) / \sqrt(t) dt, (2a)
       S2 = \sqrt(1/2/pi) \int_0^x sin(t) / \sqrt(t) dt (2b)

The integrals are calculated as follows:
- Values of X in the interval [-5,5] are looked up in Table 7.7 in [1] and interpolated ('cubic'), if necessary. This table has values of the Fresnel integrals with 7 significant digits and even linear interpolation gives an error of no more than 3e-004.
- Values outside the interval [-5,5] are evaluated using the approximations under Table 7.7 in [1]. The error is less than
3e-007.
 
NOTE: The Tables were OCR'ed and although thoroughly checked, there might be some mistakes. Please let me know if you find any.

REFERENCE
[1] Abramowitz, M. and Stegun, I. A. (Eds.). "Error Function and Fresnel Integrals." Ch. 7 in "Handbook of mathematical functions with formulas, Graphs, and mathematical tables", 9th printing. New York: Dover, pp. 295-329, 1970.

EXAMPLES OF USAGE
>> Y = fresnel(X,'c');
returns the Fresnel C integral according to eq. (0a)
>> Y = fresnel(X,'s',2);
returns the Fresnel S integral values according to eq. (2b)
>> [C,S] = fresnel(X,[],1)
returns the both the Fresnel C and S integral values according to eqs. (1a), (1b)
>> [C,S] = fresnel(X,[])
returns the both the Fresnel C and S integral values according to eqs. (0a), (0b)

• MATLAB release: MATLAB 7.6 (R2008a)

[Peter Volegov]

Errors in the S table:

0.4512445 should be 0.4572445
0.4351898 should be 0.4357898
...

Looks like 7 is often OCRed as 1. Check against fcs (File ID: 6580).

[William Sailor]

This program is a factor of ten slower than the submissions by either Barrowes or Telasula. I am not sure that reading in the tables of Abramowitz and Stegun and interpolating is a good idea, either. The results from this code differs in typically the 5th decimal place vs the other two codes, which agree well with each other.

[Christos Saragiotis]

Dear William,

1. Regarding the "slowness issue":

This function is in fact FASTER than the submissions by Telasula and Barrows (btw Barrowes' submission is not something you can run, unless you correct it).

For example for an input of 200,000 samples, in my machine I get:
      Barrows: 5.5 sec, Telasula: 3.1 sec, this function: 1.8 sec

This function is indeed slower than Barrowes' and Talasula's ones when you have inputs with very few samples (say 20 samples or so) but in that case the time of execution is of the order of msec. I think that 3 msec vs. 0.3 msec is not a waste of anyone's time.

2. Regarding the "accuracy issue":

Abramowitz's and Stegun's tables are accurate to the 10th decimal. Interpolation for other values was not my idea it was Abramowitz's and Stegun's idea. The accuracy of the interpolation was commented in the description of the file above. I don't see how your "findings" on the accuracy contribute anything new. By the way, what is the accuracy of Barrowes' or Telasula's submissions?

In my opinion it is very unfortunate that you gave a bad rating based
 - on something which doesn't hold ("slowness") and
 - on a limitation ("accuracy") which is however stated in both the description and the help of the function.

[John D'Errico]

This certainly is worth more than only the 1 star claimed by William Sailor, but it can be improved.

The help and error checking are completely satisfactory in my opinion, and an H1 line is provided.

The idea of interpolating the tables can work, but better tables, as well as a better choice of interpolant are necessary. (This code uses pchip, i.e., interp1, with the 'cubic' option, so only a C1 interpolant.) Also, it is important to note that the interpolation error will grow towards the upper end if equally spaced tables are employed.

First, I'll look at the accuracy issue.

n = 10000;
T = sort(rand(n,1)*10);
FresnelCObj = @(t) cos(pi*t.^2/2);

FCquad = zeros(size(T));
for i = 1:n
  FCquad(i) = quadgk(FresnelCObj,0,T(i),'abstol',1e-16);
end
FCpred = fresnel(T,'c',0);

max(FCquad - FCpred)
ans =
      0.000354015139247155

min(FCquad - FCpred)
ans =
     -0.000445656718834175

So in fact, the maximal error of interpolation is significant, and apparently as large as 4e-4 in places. I have verified that quadgk is indeed producing the correct results to within a precision that is close to that requested.

A plot of the interpolation error shows part of the problem, i.e., the use of the interp1 'cubic' option, which is only a C1 interpolant. The errors produced are classically what I would expect to see. One problem is that pchip is simply not designed to interpolate oscillatory curves as we have here.

plot(T,FCquad - FCpred,'.','markersize',4)
grid on
title 'Fresnel prediction errors on test data'

The author can gain an extra digit or so of accuracy in my tests simply by changing to the 'spline' option with interp1. Of course, that would tend to slow down the code slightly too, but that problem too can be resolved with careful coding. Additional accuracy can also be gained by improving the tables used for interpolation, and more accuracy yet obtained by a careful choice of how to build the interpolation tables.

This brings up the speed issue. The fresnel code is indeed reasonably fast, but the issue is whether you call it with one isolated point or a million. Since fresnel is vectorized, it is extremely efficient when called on a large set of points.

For one point on my CPU, I get a time (as predicted by timeit by Steve Eddins) of:

timeit(@() fresnel(pi,'c',0))
ans =
            0.000694679339

In this test, T is a vector with 1e6 elements in it.

timeit(@() fresnel(T,'c',0))
ans =
            0.266534082114

.26653/.000686
ans =
          388.527696793003

Thus the time required per point on the large sample is now only roughly 0.00000027 seconds per point. Well vectorized code pays dividends.

What should I rate this? Since the accuracy can be improved with little effort on the part of the author, I'll give it 4 stars.

[Christos Saragiotis]

John,

thank you very much for your comments and the time you spent on thoroughly checking this function.

I will try to incorporate your suggestions, when I find some time