function [fs,fc,ff,fg]=fresnel(x,xc)
%FRESNEL(x) Calculate Fresnel sin and cos integrals
% x is the argument to the fresnel function
% xc (default=1.8) is the cutoff which determines the method used
%                  to calculate the integral
%
% for |x|>xc ff and fg are the Auxiliary functions
% for |x|>xc ff and fg = NaN
%
% if fresnel is called with no input arguments a self test is run
%
%       Fresnel Integrals
%
%         /x
% fc(x) = |   cos(pi/2*x^2) dx                                 (7.3.1)
%         /0
%
%         /x
% fs(x) = |   sin(pi/2*x^2) dx                                 (7.3.2)
%         /0
%
%    Auxiliary functions
%
%    ff = (1/2-fs) * cos(pi/2x^2) - (1/2-fc) * sin(pi/2x^2)    (7.3.5)
%
%    fg = (1/2-fc) * cos(pi/2x^2) + (1/2-fs) * sin(pi/2x^2)    (7.3.6)
%
%  d fs
%  ---- =  sin(pi/2*x^2)
%  d x
%
%
%  d fc
%  ---- =  cos(pi/2*x^2)
%  d x
%
%  for |x| < xc uses a powers series about x=0         
%               accuracy is better than 5e-16
%  for |x| > xc uses a minimax rational approximation
%               based on the auxilliary functions f and g
%               accuracy is better than 1e-9
%
% The accuracy can be checked by numerically calculating
% dfs/dx and dfc/dx. For all x this is better than 2e-8 
% 
% To test the accuracy call with no arguments
% Approximations and 20 digit accurate results were generated with maple 6
%
% Definitions for the Fresnel integrals and auxialliary functions
% f and g are taken from 
%      Handbook of Mathematical Functions
%      Abramowitz and Stegun (9th printing 1970)
%      Section 7.3 p300
%
% uses power series for x<xc=1.8                  accurate to machine precision
% uses minimax rational approximation for x>xc=1.8     accurate to about 1e-9
%
% 
% Released under the GPL
% jnm11@amtp.cam.ac.uk
% J. N. McElwaine
% Department of Applied Mathematcis and Theoretical Physics
% University of Cambridge 
% Silver Street
% UK
% Versions
%
% 1.0 October 2001
% 
% 1.01 29th November 2001
%
%
% TODO
% The accuracy should be improved to machine precision for larger values of x
%

if nargin==0
  test_fresnel;
  return;
end

% Cut point for switching between power series and rational approximation
if nargin<2
  xc=[];
end
if isempty(xc)
  xc=1.8; 
end

sgnx = sign(x);
x = abs(x);

f1 = find(x <= xc);
f2 = find(x >  xc);

s = pi/2*x(f1).^2;
ss = s.^2;
is = sqrt(2/pi)./x(f2);

fs = repmat(NaN,size(x));
fc = fs;
ff = fs;
fg = fs;

% Approximations for x < 1.6 accurate to eps

fs(f1)=x(f1).*s.*(1/3+(-.23809523809523809524e-1+(.75757575757575757576e-3+(-.13227513227513227513e-4+(.14503852223150468764e-6+(-.10892221037148573380e-8+(.59477940136376350368e-11+(-.24668270102644569277e-13+(.80327350124157736091e-16+(-.21078551914421358249e-18+(.45518467589282002862e-21+(-.82301492992142213568e-24+(.12641078988989163522e-26+(-.16697617934173720270e-29+.19169428621097825308e-32.*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss); 

fc(f1)=x(f1).*(1+(-.1+(.46296296296296296296e-2+(-.10683760683760683761e-3+(.14589169000933706816e-5+(-.13122532963802805073e-7+(.83507027951472395917e-10+(-.39554295164585257634e-12+(.14483264643598137265e-14+(-.42214072888070882330e-17+(.10025164934907719167e-19+(-.19770647538779051748e-22+(.32892603491757517328e-25+(-.46784835155184857737e-28+.57541916439821717722e-31.*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss).*ss);

ff(f2)=sqrt(2/pi)*is.*(.36634901025764670680+(.84695666222194589182+(-3.7301273487349839902+(10.520237918558456228+(-10.472617402497801126+4.4169634834911107719.*is).*is).*is).*is).*is)./(.73269802661202980832+(1.6939102288852613619+(-7.4599994789665215344+(21.032436583848862358+(-20.269535575486282590+8.9995024877628789836.*is).*is).*is).*is).*is);

fg(f2)=sqrt(2/pi)*is.^3.*(.10935692320079194659+(.72025533055541994934+(-2.5630322339412334317+(7.2404268469720856874+(-7.0473933910697823445+2.9315207450903789503.*is).*is).*is).*is).*is)./(.43742772185339366619+(2.8810049959848098312+(-10.250672312139277005+(28.912791022417600679+(-25.740131167525284201+15.145134363709932380.*is).*is).*is).*is).*is);

s = pi/2*x(f2).^2;
cx = cos(s);
sx = sin(s);
fs(f2) = 1/2 - ff(f2).*cx - fg(f2).*sx;
fc(f2) = 1/2 + ff(f2).*sx - fg(f2).*cx;

fs = fs.*sgnx;
fc = fc.*sgnx;


function test_fresnel

x1 = 0:0.1:1.6;
s1 = [0
      .52358954761221059949e-3
      .41876091616567616310e-2
      .14116998006576585807e-1
      .33359432660613180394e-1
      .64732432859999277611e-1
      .11054020735938696133
      .17213645786347745336
      .24934139305391778373
      .33977634439314021465
      .43825914739035476608
      .53649791109682043483
      .62340091854624967227
      .68633328553465011378
      .71352507736341211296
      .69750496008209301308
      .63888768350938090306]';

c1 = [0
      .99997532627085068050e-1
      .19992105759445308522
      .29940097605204721038
      .39748075917235943981
      .49234422587144639288
      .58109544699165232707
      .65965235190451039092
      .72284417189635611829
      .76482302127332649986
      .77989340037682282947
      .76380666606201199073
      .71543772292307339595
      .63855045472702925725
      .54309578354625638856
      .44526117603982153506
      .36546168344048770958]';

x2=[1.6,1.7,1.8,1.9,2,2.2,2.4,2.6,2.8,3.0,3.5,4.5,5,6,7,8,9,10,15,20,50];

% map(Fresnelf[1.61.71.81.92.2.22.42.62.83.03.54.55.6.7.8.9.10.15.20.50.]);
f1=[.19219389591936835161 
    .18200800122326605067
    .17274483591298635463 
    .16430085282185039519
    .15658432163630175780 
    .14302018320281828056
    .13151755794910432090 
    .12166560172918624682
    .11314822051641228878 
    .10572078929768562956
    .090765558315310835015 
    .070683539588491872607
    .063631188704012231102 
    .053039238763069722376
    .045467092546969810327 
    .039785785606985516138
    .035366127468119852551 
    .031830021415117759597
    .021220531674373457935 
    .015915464074046107285
    .0063661974140612585467]';

% map(Fresnelg[1.61.71.81.92.2.22.42.62.83.03.54.55.6.7.8.9.10.15.20.50.]);

g1=  [.021256848886994052896
      .018174092917668532180
      .015627825544284288857
      .013512944617910355890
      .011746593924659245500
      .0090108056825678170591
      .0070410079901931512250
      .0055940626784959684812
      .0045112596246994357667
      .0036870010326249639024
      .0023401756317289203885
      .0011078332746661747626
      .00080861808288311324807
      .00046853214449887981172
      .00029521054655322413172
      .00019781962280286444302
      .00013895437031507013829
      .00010130579448427638585
      .000030020190297256575884
      .000012665027655611812994
      .81056927203204418984e-6]';

%map(FresnelS,[1.6,1.7,1.8,1.9,2,2.2,2.4,2.6,2.8,3.0,3.5,4.5,5,6,7,8,9,10,15,20,50]);

s2 = [.63888768350938090306
      .54919594032156850105
      .45093876926758310142
      .37334731781698114416
      .34341567836369824220
      .45570461212465689284
      .61968996494568358267
      .54998932315271946901
      .39152844354317181852
      .49631299896737503610
      .41524801197243752390
      .43427297504870358947
      .49919138191711688675
      .44696076123693027762
      .49970478945344677587
      .46021421439301448386
      .49986104562968492986
      .46816997858488224040
      .49996997980970274342
      .48408453592595389271
      .49363380258593874145]';

%map(FresnelC,[1.6,1.7,1.8,1.9,2,2.2,2.4,2.6,2.8,3.0,3.5,4.5,5,6,7,8,9,10,15,20,50]);

c2 = [.36546168344048770958
      .32382687600390025374
      .33363292722155710067
      .39447053489152294822
      .48825340607534075450
      .63628604490331945781
      .55496140585642812795
      .38893749619196902786
      .46749165169890597878
      .60572078929768562956
      .53257243502800084533
      .52602591505353874137
      .56363118870401223110
      .49953146785550112019
      .54546709254696981033
      .49980218037719713556
      .53536612746811985255
      .49989869420551572361
      .52122053167437345793
      .49998733497234438819
      .49999918943072796796]';

[fs fc ff fg] =  fresnel(x1);

subplot(3,1,1)
plot(x1,fs-s1,x1,fc-c1);
title('Fresnel integrals errors for small x');
subplot(3,1,2)
[fs fc ff fg] =  fresnel(x2);
plot(x2,fs-s2,x2,fc-c2);
title('Fresnel integrals errors for large x');

subplot(3,1,3)
x=linspace(0,10,1000);
dx = eps+x*sqrt(eps);
[sp cp] =  fresnel(x+dx);
[sm cm] =  fresnel(x-dx);
plot(x,cos(pi/2*x.^2)-(cp-cm)./(2*dx));
plot(x,sin(pi/2*x.^2)-(sp-sm)./(2*dx));
title('Differential of Fresnel integrals errors for all x');


return