File Exchange

image thumbnail

FresnelS and FresnelC

version 1.1.0.0 (2.34 MB) by John D'Errico
Efficient and accurate computation of the Fresnel sine and cosine integrals

8 Downloads

Updated 03 May 2012

View Version History

View License

I noticed the many codes on the FEX to compute the Fresnel integrals for real arguments, and it left me wondering how I might try solving this problem in MATLAB for both high accuracy and high efficiency.

The approach I took yields a maximum error of roughly 1e-14 as far as I could get reasonable values to compare it to. (The screenshot shows the predicted error for a sampling of points.)

I've supplied functions for both the Fresnel sine and cosine integrals, as well as a .pdf file that explains the approach I took.

Evaluate the Fresnel cosine integral C(x) at x = 1.38

>> fresnelC(1.38,0)
ans =
0.562975925772444

Verify the correctness of this value using quadgk.

>> FresnelCObj = @(t) cos(pi*t.^2/2);
>> quadgk(FresnelCObj,0,1.38,'abstol',1e-15')
ans =
0.562975925772444

Now, how fast is fresnelC? Using Steve Eddins timeit code to yield an accurate estimate of the time required, we see that it is reasonably fast for scalar input.

>> timeit(@() fresnelC(1.38))
ans =
0.000193604455833333

More importantly, these functions are properly vectorized. So 1 million evaluations are easy to do, and are much faster than 1 million times the time taken for one evaluation.

>> T = rand(1000000,1);
>> tic
>> FCpred = fresnelC(T);
>> toc
Elapsed time is 0.226884 seconds.

Cite As

John D'Errico (2020). FresnelS and FresnelC (https://www.mathworks.com/matlabcentral/fileexchange/28765-fresnels-and-fresnelc), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (8)

Patrick Hew

Sorry, just a quick edit to my previous comment. The correct footprints should be
FCint = fresnelC0(X)
FSint = fresnelS0(X)
to calculate the type 0 integrals where X is non-negative. It is for the caller to keep track of which entries in X should be flipped back to negative (or not).

Patrick Hew

This submission works very nicely, thank you for providing it.

As a suggestion for a future upgrade (if any), perhaps the "core" code from lines 112 onwards could be refactored into functions:
FCint = fresnelC0(XS,S)
FSint = fresnelS0(XS,S)
where fresnelC0 & fresnelS0 calculate the type 0 integrals for X where X(k) = XS(k) if S(k) >= 0 and X(k) = -XS(k) if S(k) < 0. The original fresnelC & fresnelS functions would delegate to fresnelC0 & fresnelS0, but users would also have the option of calling fresnelC0 & fresnelS0 directly.

The reason is that in profiling, something like 35 percent of the time is spent on identifying the type of integral (line 85, checking fresnelType), and another 10 percent in handling the polarity of X (line 102 to flip X(S) to positive and then line 154 to flip FCint(S) to negative).

Andreas Trzuskowsky

Great work!
I needed it in a Simulink Model within a „Matlab Function“ block. So I had to modify it a bit in order to make it ready for code generation. Unfortunately Simulink states “PP forms constructed in MATLAB cannot be used directly by PPVAL for code generation. Use UNMKPP to extract the data and MKPP to create a MATLAB Coder compatible PP form.” Therefore, after loading the PP-data using the coder.load command, the data needs to be split and conducted again.
So I replaced line 131

load _Fresnel_data_ FCspl

with

content = coder.load('_Fresnel_data_.mat');
[breaks,coefs,~,~,~] = unmkpp(content.FCspl);
FCspl = mkpp(breaks, coefs);

Accordingly for FSspl. In addition, since code generation does not support warnings the line

coder.extrinsic('warning');

needs to be added to make the functions usable for code generation.

gzrollins

Very useful functions. Very well documented. First class work! Thank you!

Amir Khabbazi Oskouei

John D'Errico

New version submitted - thanks to Felipe.

Felipe

Hi John. Thanks for citing ("acknowledging") related submissions. There are two new ones, that came after yours: 33577 and 34134. You might want to cite these, too. That way folks will find your submission in all cases. I'm trying to kill duplicates. I'm assuming yours is superior -- both in accuracy and speed, not to mention readability.

John Kot

MATLAB Release Compatibility
Created with R2010a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!