'--------------------------------------------------------------- ' TESTFRES.BAS ' ' http://www.achilles.net/~jtalbot/amateur/diffract ' ' by John Talbot, jtalbot@achilles.net Apr 19,1997. ' ' This version works on any computer platform ' ' Test the accuracy and speed of various methods to compute Fresnel ' integrals against tabulated values of the fresnel integrals from ' p.321 in HANDBOOK OF MATHEMATICAL FUNCTIONS by Abramowitz and Stegun ' (Fresnel3 was chosen to be the most accurate method) ' '--------------------------------------------------------------- DEFDBL A-Z DIM SHARED A,B,PI,B1(11),B2(11),B3(11),B4(11),j(100) PI = 3.141592653589# '--------------------------------------------------------------- TestFresnel: N=16 'GOSUB Terms GOSUB Coefficients RESTORE FresnelData FOR i=1 TO 1000 READ X,C0,S0 :IF X=999 THEN END U=X*X*PI/2 GOTO Skip IF X<4 THEN FresnelSmall C1,S1,U ELSE FresnelLarge C1,S1,X END IF Fresnel2 C2,S2,X GetFresnel C3,S3,U,N Skip: Fresnel3 C4,S4,X PRINT USING "##.####### ";C0,S0; PRINT USING " X= ##.### U= ####.### ";X,U 'PRINT USING "##.####### ";C1,S1 'PRINT USING "##.####### ";C2,S2 'PRINT USING "##.####### ";C3,S3 PRINT USING "##.####### ";C4,S4 PRINT USING "##.####### ";C0-C4,S0-S4 'IF C0>0 AND S0>0 THEN ' F=3300 'pixels on a laser printer ' PRINT USING " #####";F*(C4-C0)/C0,F*(S4-S0)/S0 'END IF PRINT NEXT i RETURN SUB fresnel(C,S,XOrigin,XPos) STATIC U = A*(B*XOrigin-XPos) IF U<4 THEN FresnelSmall C,S,U*U*PI/2 ELSE FresnelLarge C,S,U END IF 'IF U<0 THEN C=-C : S=-S END SUB '----------------------------------------------------- ' Asymptotic expansion of Fresnel integral for large X ' see Abramowitz p.301 ' Accurate for 4 < X < Infinity SUB FresnelLarge(C,S,X) STATIC 'will this work for X<0 ??? Try it out ... X2=X*X F = ( 1# + .926#*X ) / ( 2# + 1.792#*X + 3.104#*X2 ) G = 1# / ( 2# + 4.142#*X + 3.492#*X2 + 6.67#*X2*X ) U = PI*X2/2 CosU = COS(U) SinU = SIN(U) C = .5# + F*SinU - G*CosU S = .5# - F*CosU - G*SinU END SUB '------------------------------------------------- 'SUBROUTINE: Solve The Fresnel Integrals S2(x),C2(x) ' stable for 0 < x < 25 , SQR(2*x/PI) < 4.1 ' can use up to 49 terms without convergeance problems SUB FresnelSmall(C,S,X) STATIC IF X=0 THEN S=0:C=0:EXIT SUB '32 IF X>32 THEN 'Aproximation: Max Deviation plus or minus 2E-4 C=.5+SIN(X)/SQR(2*PI*X) S=.5-COS(X)/SQR(2*PI*X) EXIT SUB END IF Sign=1:PowerX=1:B=1: Factorial=1:Fac=1 C=0:S=0:LastC=0:LastS=0 FOR k=1 TO 85 'do not change 85 (Max before overflow) dc=Sign*PowerX/(B*Factorial) : C=C+dc Factorial=Fac*Factorial : Fac=Fac+1 PowerX=X*PowerX B=B+2 ds=Sign*PowerX/(B*Factorial) : S=S+ds Factorial=Fac*Factorial : Fac=Fac+1 PowerX=X*PowerX B=B+2 Sign=-Sign IF LastS=S AND LastC=C THEN NoMoreTerms 'for a better control of tolerance:use ds,dc LastS=S:LastC=C NEXT k NoMoreTerms: B=SQR(X*2/PI) C=C*B S=S*B END SUB '----------------------------------------------------------------- 'Numerical evaluation of Fresnel integrals using taylor series '(from the original FORTRAN program, doesn't seem to work yet ! ) SUB Fresnel2(C,S,X) STATIC U = X*X*PI/2 IF U>4 THEN XF=4/U : C=0 : S=0 : Z=1 FOR i=0 TO 8 C = C + B3(i)*Z S = S + B4(i)*Z Z = Z*XF NEXT i Z=SQR(U) : D=COS(Z) : S=SIN(Z) C = .5 + Z*(D*A + S*B) S = .5 + Z*(S*A - D*B) ELSE XF=U : C=0 : S=0 : Z=1 FOR i=0 TO 6 C = C + B1(i)*Z S = S + B2(i)*Z Z = Z*XF NEXT i Z=SQR(U) : C=Z*C : S=X*Z*S END IF IF X<0 THEN C =-C : S =-S END SUB '---------------------------------------------------------------- 'Read coefficients of the power series for evaluating the fresnel 'intergrals. Terms: DATA 1.840962E-1 , 8.02649E-2 , 1.994712E-1 , 4.444091E-9 DATA 1.102544E-2 , 6.1213E-3 , -1.207998E-6 , -2.493322E-2 DATA 1.020418E-3 , 2.441816E-4 , -9.314911E-3 , -1.606428E-5 DATA 3.273308E-5 , 5.051141E-6 , -4.027145E-4 , 5.972151E-3 DATA 5.451182E-7 , 5.883158E-8 , 7.428246E-3 , -3.095341E-4 DATA 5.244297E-9 , 6.677681E-10 , -7.27169E-2 , -6.792801E-3 DATA 5.100785E-11 , 0 , 3.401409E-3 , 7.970943E-3 DATA 0 , 0 , -6.633926E-4 , -4.169289E-3 DATA 0 , 0 , 0 , 8.768258E-4 RESTORE Terms FOR i=0 TO 8 READ B1(i),B2(i),B3(i),B4(i) NEXT i RETURN ' tabulated values of the fresnel integrals from p.321 in ' HANDBOOK OF MATHEMATICAL FUNCTIONS by Abramowitz and Stegun FresnelData: DATA 0.0, 0.0000000, 0.0000000 DATA 0.1, 0.0999975, 0.0005236 DATA 0.2, 0.1999211, 0.0041876 DATA 0.3, 0.2994010, 0.0141170 DATA 0.4, 0.3974808, 0.0333594 DATA 0.5, 0.4923442, 0.0647324 DATA 0.6, 0.5810954, 0.1105402 DATA 0.7, 0.6596524, 0.1721365 DATA 0.8, 0.7228442, 0.2493414 DATA 0.9, 0.7648230, 0.3397763 DATA 1.0, 0.7798934, 0.4382591 DATA 1.1, 0.7638067, 0.5364979 DATA 1.2, 0.7154377, 0.6234009 DATA 1.3, 0.6385505, 0.6863333 DATA 1.4, 0.5430958, 0.7135251 DATA 1.5, 0.4452612, 0.6975050 DATA 1.6, 0.3654617, 0.6388877 DATA 1.7, 0.3238269, 0.5491960 DATA 1.8, 0.3336329, 0.4509388 DATA 1.9, 0.3944705, 0.3733473 DATA 2.0, 0.4882534, 0.3434157 DATA 2.1, 0.5815641, 0.3742734 DATA 2.2, 0.6362860, 0.4557046 DATA 2.3, 0.6265617, 0.5531516 DATA 2.4, 0.5549614, 0.6196900 DATA 2.5, 0.4574130, 0.6191818 DATA 2.6, 0.3889375, 0.5499893 DATA 2.7, 0.3924940, 0.4529175 DATA 2.8, 0.4674917, 0.3915284 DATA 2.9, 0.5623764, 0.4101406 DATA 3.0, 0.6057208, 0.4963130 DATA 3.1, 0.5615939, 0.5818159 DATA 3.2, 0.4663203, 0.5933495 DATA 3.3, 0.4056944, 0.5192861 DATA 3.4, 0.4384917, 0.4296495 DATA 3.5, 0.5325724, 0.4152480 DATA 3.6, 0.5879533, 0.4923095 DATA 3.7, 0.5419457, 0.5749804 DATA 3.8, 0.4480949, 0.5656187 DATA 3.9, 0.4223327, 0.4752024 DATA 4.0, 0.4984260, 0.4205158 DATA 4.1, 0.5736956, 0.4757983 DATA 4.2, 0.5417192, 0.5631989 DATA 4.3, 0.4494412, 0.5539959 DATA 4.4, 0.4383329, 0.4622680 DATA 4.5, 0.5260259, 0.4342730 DATA 4.6, 0.5672367, 0.5161923 DATA 4.7, 0.4914265, 0.5671455 DATA 4.8, 0.4337966, 0.4967502 DATA 4.9, 0.5001610, 0.4350674 DATA 5.0, 0.5636312, 0.4991914 DATA 999, 0 , 0 '------------------------------------------------------- SUB GetFresnel(C,S,X,N) STATIC N%=N 'Starting order for backward recursion N1%=.8*N 'Upper bessel order to start Fresnel Sums j(N%) =0 j(N%-1)=1E-12 C=0 : S=0 R2X=2/X : RX=1/X FOR i%=N%-1 TO 1 STEP -1 j(i%-1)=(R2X*i%+RX)*j(i%)-j(i%+1) PRINT "N=";i%; IF i%