function [C,S] = fresnel(x,cs,var)
% 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,
% S = \int_0^x sin(pi*t^2/2) dt
% Also, it evaluates the following variations of the Fresnel integrals
% C1 = \sqrt(2/pi) \int_0^x cos(t^2) dt,
% S1 = \sqrt(2/pi) \int_0^x sin(t^2) dt
% and
% C2 = \sqrt(1/2/pi) \int_0^x cos(t) / \sqrt(t) dt,
% S2 = \sqrt(1/2/pi) \int_0^x sin(t) / \sqrt(t) dt
%
% SYNTAX
% Y = fresnel(X,'TYPE',VAR);
% Y = fresnel(X,'TYPE');
% [C,S] = fresnel(X,[],VAR)
% [C,S] = fresnel(X,[])
% where,
% X is an array of values on which the Fresnel integral are to
% be evaluated,
% 'TYPE' is either 'c' or 's' for Fresnel C integral or Fresnel S
% integral. Default is [], which returns both Fresnel C and S.
% VAR is either 0,1 or 2 and refers to the variations described
% above (default is 0)
%
% The output(s) have the same shape as the input.
%
% HOW THE INTEGRALS ARE CALCULATED
% - 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, so there might be some mistakes.
%
% REFERENCES
% [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.
%
% By C. Saragiotis, May 2008
%
%% Initialization
if nargin < 3
var = 0;
if nargin < 2
cs='both';
if nargin<1
error('Function FRESNEL must have at least one argument');
end
end
end
if isempty(var), var=0; end
%% Inputs check
s = size(x);
x = reshape(x,[prod(s),1]);
if var~=0 && var~=1 && var~=2,
error('Input VAR must be either 0, 1 or 2');
end
if isempty(cs), cs = 'both'; end
cs = lower(cs);
if ~strcmp(cs,'c') && ~strcmp(cs,'s') && ~strcmp(cs,'both'),
error('Input TYPE must be either ''c'' or ''s'' or []');
end
row = size(x,2)>1;
if row, x = x.'; end
if imag(x)~=0
warning('Only the real part of x is considered') %#ok<WNTAG>
x = real(x);
end
%% Main
sgn = sign(x);
x = abs(x);
if var==1,
x = sqrt(2/pi)*x;
elseif var==2,
x = sqrt(2*x/pi);
end
% Discriminate x's <= 5 and x's >5
i05 = x<=5;
x05 = x(i05);
switch cs
case 'c'
C = evalC(x,i05,x05,sgn,s);
S = [];
case 's'
C = evalS(x,i05,x05,sgn,s);
S = [];
otherwise
C = evalC(x,i05,x05,sgn,s);
S = evalS(x,i05,x05,sgn,s);
end
end
%% #### Local functions ####
%% evalC
function C = evalC(x,i05,x05,sgn,s)
C = 0.5 + (0.3183099-0.0968./x.^4) .* sin(pi/2*x.^2)./x - (0.10132-0.154./x.^4) .* cos(pi/2*x.^2)./x.^3;
C(i05) = tableC(x05); % Use table for values in [-2,2]
C = sgn.*C; % Fresnel C is an odd function
C = reshape(C,s); % Reshape to original shape
end
%% evalS
function S = evalS(x,i05,x05,sgn,s)
S = 0.5 - (0.3183099-0.0968./x.^4) .* cos(pi/2*x.^2)./x - (0.10132-0.154./x.^4) .* sin(pi/2*x.^2)./x.^3;
S(i05) = tableS(x05); % Use table for values in [-2,2]
S = sgn.*S; % Fresnel S is an odd function
S = reshape(S,s); % Reshape to original shape
end
%% tableC
function C = tableC(x)
fC05 = [
0.0000000; % x=0
0.0200000;0.0400000;0.0599998;0.0799992;0.0999975;0.1199939;0.1399867;
0.1599741;0.1799534;0.1999211;0.2198729;0.2398036;0.2597070;0.2795756;
0.2994010;0.3191731;0.3388806;0.3585109;0.3780496;0.3974808;0.4167868;
0.4359482;0.4549440;0.4737510;0.4923442;0.5106969;0.5287801;0.5465630;
0.5640131;0.5810954;0.5977737;0.6140094;0.6297625;0.6449912;0.6596524;
0.6737012;0.6870920;0.6997779;0.7117113;0.7228442;0.7331283;0.7425154;
0.7509579;0.7584090;0.7648230;0.7701563;0.7743672;0.7774168;0.7792695;
0.7798934; % x=1
0.7792611;0.7773501;0.7741434;0.7696303;0.7638067;0.7566760;0.7482494;
0.7385468;0.7275968;0.7154377;0.7021176;0.6876947;0.6722378;0.6558263;
0.6385505;0.6205111;0.6018195;0.5825973;0.5629759;0.5430958;0.5231058;
0.5031623;0.4834280;0.4640705;0.4452612;0.4271732;0.4099799;0.3938529;
0.3789596;0.3654617;0.3535120;0.3432529;0.3348132;0.3283061;0.3238269;
0.3214502;0.3212283;0.3231887;0.3273325;0.3336329;0.3420339;0.3524496;
0.3647635;0.3788293;0.3944705;0.4114824;0.4296333;0.4486669;0.4683056;
0.4882534; % x=2
0.5082004;0.5278273;0.5468106;0.5648279;0.5815641;0.5967175;0.6100060;
0.6211732;0.6299953;0.6362860;0.6399031;0.6407525;0.6387928;0.6340383;
0.6265617;0.6164945;0.6040269;0.5894065;0.5729344;0.5549614;0.5358811;
0.5161229;0.4961428;0.4764135;0.4574130;0.4396132;0.4234672;0.4093965;
0.3977791;0.3889375;0.3831273;0.3805280;0.3812350;0.3852532;0.3924940;
0.4027739;0.4158168;0.4312585;0.4486546;0.4674917;0.4872004;0.5071721;
0.5267766;0.5453821;0.5623764;0.5771878;0.5893060;0.5983019;0.6038456;
0.6057208; % x=3
0.6038373;0.5982378;0.5891011;0.5767401;0.5615939;0.5442158;0.5252553;
0.5054356;0.4855276;0.4663203;0.4485896;0.4330655;0.4204005;0.4111397;
0.4056944;0.4043199;0.4070996;0.4139366;0.4245518;0.4384917;0.4551437;
0.4737596;0.4934870;0.5134062;0.5325724;0.5500611;0.5650132;0.5766802;
0.5844643;0.5879533;0.5869464;0.5814710;0.5717875;0.5583818;0.5419457;
0.5233449;0.5035770;0.4837194;0.4648719;0.4480949;0.4343486;0.4244343;
0.4189443;0.4182216;0.4223327;0.4310568;0.4438917;0.4600770;0.4786351;
0.4984260; % x=4
0.5182154;0.5367505;0.5528404;0.5654347;0.5736956;0.5770588;0.5752776;
0.5684474;0.5570075;0.5417192;0.5236206;0.5039608;0.4841163;0.4654961;
0.4494412;0.4371250;0.4294640;0.4270439;0.4300679;0.4383329;0.4512359;
0.4678105;0.4867941;0.5067195;0.5260259;0.5431811;0.5568046;0.5657827;
0.5693657;0.5672367;0.5595481;0.5469186;0.5303913;0.5113538;0.4914265;
0.4723271;0.4557230;0.4430830;0.4355428;0.4337966;0.4380247;0.4478669;
0.4624440;0.4804290;0.5001610;0.5197951;0.5374734;0.5515025;0.5605194;
0.5636312]; % x=5
C = interp1(0:0.02:5, fC05,x,'cubic');
end
%% tableS
function S = tableS(x)
fS05 = [
0.0000000; % x=0
0.0000042;0.0000335;0.0001131;0.0002681;0.0005236;0.0009047;0.0014367;
0.0021444;0.0030531;0.0041876;0.0055730;0.0072340;0.0091954;0.0114816;
0.0141170;0.0171256;0.0205311;0.0243568;0.0286255;0.0333594;0.0385802;
0.0443085;0.0505642;0.0573663;0.0647324;0.0726789;0.0812206;0.0903708;
0.1001409;0.1105402;0.1215759;0.1332528;0.1455729;0.1585354;0.1721365;
0.1863689;0.2012221;0.2166816;0.2327288;0.2493414;0.2664922;0.2841498;
0.3022780;0.3208355;0.3397763;0.3590493;0.3785981;0.3983612;0.4182721;
0.4382591; % x=1
0.4582458;0.4781508;0.4978884;0.5173686;0.5364979;0.5551792;0.5733128;
0.5907966;0.6075274;0.6234009;0.6383134;0.6521619;0.6648456;0.6762672;
0.6863333;0.6949562;0.7020550;0.7075567;0.7113977;0.7135251;0.7138977;
0.7124878;0.7092816;0.7042812;0.6975050;0.6889888;0.6787867;0.6669713;
0.6536346;0.6388877;0.6228607;0.6057026;0.5875804;0.5686783;0.5491960;
0.5293473;0.5093584;0.4894649;0.4699094;0.4509388;0.4328006;0.4157397;
0.3999944;0.3857925;0.3733473;0.3628537;0.3544837;0.3483830;0.3446665;
0.3434157; % x=2
0.3446748;0.3484487;0.3547004;0.3633498;0.3742734;0.3873037;0.4022309;
0.4188045;0.4367363;0.4557046;0.4753585;0.4953241;0.5152111;0.5346203;
0.5531516;0.5704128;0.5860284;0.5996489;0.6109596;0.6196900;0.6256211;
0.6285938;0.6285143;0.6253598;0.6191818;0.6101076;0.5983406;0.5841575;
0.5679042;0.5499893;0.5308753;0.5110679;0.4911035;0.4715352;0.4529175;
0.4357898;0.4206603;0.4079890;0.3981724;0.3915284;0.3882841;0.3885643;
0.3923850;0.3996480;0.4101406;0.4235387;0.4394139;0.4572445;0.4764306;
0.4963130; % x=3
0.5161942;0.5353629;0.5531195;0.5688028;0.5818159;0.5916511;0.5979129;
0.6003366;0.5988034;0.5933495;0.5841697;0.5716147;0.5561806;0.5384935;
0.5192861;0.4993695;0.4796004;0.4608446;0.4439382;0.4296495;0.4186411;
0.4114369;0.4083928;0.4096754;0.4152480;0.4248672;0.4380883;0.4542817;
0.4726592;0.4923095;0.5122412;0.5314321;0.5488815;0.5636638;0.5749804;
0.5822056;0.5849261;0.5829692;0.5764191;0.5656187;0.5511574;0.5338432;
0.5146622;0.4947245;0.4752024;0.4572613;0.4419892;0.4303279;0.4230117;
0.4205158; % x=4
0.4230199;0.4303900;0.4421781;0.4576445;0.4757983;0.4954571;0.5153214;
0.5340587;0.5503941;0.5631989;0.5715723;0.5749103;0.5729547;0.5658205;
0.5539959;0.5383155;0.5199077;0.5001173;0.4804108;0.4622680;0.4470706;
0.4359933;0.4299086;0.4293116;0.4342730;0.4444234;0.4589736;0.4767689;
0.4963756;0.5161923;0.5345797;0.5499967;0.5611328;0.5670244;0.5671455;
0.5614619;0.5504452;0.5350416;0.5165982;0.4967502;0.4772800;0.4599575;
0.4463774;0.4378082;0.4350674;0.4384348;0.4476156;0.4617567;0.4795178;
0.4991914]; % x=5
S = interp1(0:0.02:5, fS05,x,'cubic');
end
