1 files changed, 74 insertions, 0 deletions
diff --git a/transverse/AiryE.m b/transverse/AiryE.m
new file mode 100644
index 0000000..b03b5ed
--- /dev/null
+++ b/transverse/AiryE.m
@@ -0,0 +1,74 @@
+%---------------------------------------------------------------
+% PROGRAM: AiryE
+% AUTHOR:  Andri M. Gretarsson
+% DATE:    8/7/03
+%
+%
+% SYNTAX: E=AiryE([a,w0,L,lambda <,N>],r);
+%              <,...> indicates optional argument
+%
+% Returns the Field of an Airy disk as a function of radius.  The waist
+% of the beam is assumed to be at the aperture.
+%
+% a      = radius of aperture
+% w0     = radius of Gaussian beam at aperture
+% L      = propagation distance (axially) from aperture
+% lambda = wavelength
+% N      = order to which to carry the calculation
+%
+% Last updated: 8/7/03 by AMG
+%
+%---------------------------------------------------------------
+%% SYNTAX: E=AiryE([a,w0,L,lambda <,N>],r);
+%---------------------------------------------------------------
+
+function E=AiryE(params,r)
+
+a=params(1);
+w0=params(2);
+L=params(3);
+lambda=params(4);
+if length(params)>=5, Norders=params(5); else Norders=2; end
+if Norders<1, Norders=1; end
+
+%N=a^2/L/lambda;
+%q0=i*pi*w0.^2/lambda;
+%q=q0+L;
+%E=q0/q .* exp(-i*pi*N*(r./a).^2) .* ( 1 - exp(-i*pi*N-a^2/w0^2) .* besselj(0,2*pi*N*r./a) );
+
+k=2*pi/lambda;
+F=(a/w0)^2/pi;
+R=sqrt(r.^2+L^2);
+eta=pi*F/a^2-i*k/2./R;
+gamma=-k.*r./R;
+
+
+Sum1=0; 
+Sum2=0;
+for m=0:Norders-1
+    Sum1=Sum1+(-gamma/2/a./eta).^m .* besselj(m,a.*gamma);
+    Sum2=Sum2+(gamma/2/a./eta).^m + besselj(1-m,a*gamma);
+end
+
+
+A = -i*k*L/2./eta./R.^2 .* exp(-i*k.*R);
+B = exp(-gamma.^2/4./eta)  -  exp(-eta*a^2) .* Sum1;
+C = gamma/a./eta.*exp(-gamma.^2/4./eta)  -  2*exp(-eta*a^2) .* Sum2;
+
+E= sqrt(2)/pi/w0 * A .* (B + (i*a*r./R.^2).*C);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+