diff options
| author | Eugeniy Mikhailov <evgmik@gmail.com> | 2014-05-13 15:35:26 -0400 |
|---|---|---|
| committer | Eugeniy Mikhailov <evgmik@gmail.com> | 2014-05-13 15:58:56 -0400 |
| commit | d1b62466d81b164656588ac1c55acff15ee9ea43 (patch) | |
| tree | cce78f90b0768361c4a268b946848842e1e47456 /transverse | |
| download | optics_toolkit-d1b62466d81b164656588ac1c55acff15ee9ea43.tar.gz optics_toolkit-d1b62466d81b164656588ac1c55acff15ee9ea43.zip | |
initial
The optics_toolkit code taken from
http://mercury.pr.erau.edu/~greta9a1/downloads/index.html
the older version is also available at mathwork web site
http://www.mathworks.com/matlabcentral/fileexchange/15459-basic-paraxial-optics-toolkit
Diffstat (limited to 'transverse')
| -rw-r--r-- | transverse/AiryE.m | 74 | ||||
| -rw-r--r-- | transverse/AiryI.m | 25 | ||||
| -rw-r--r-- | transverse/AiryI_Fraunhofer.m | 33 | ||||
| -rw-r--r-- | transverse/HermiteGaussianE.m | 134 | ||||
| -rw-r--r-- | transverse/LaguerreGaussianE.m | 189 | ||||
| -rw-r--r-- | transverse/LaguerrePoly.m | 24 | ||||
| -rw-r--r-- | transverse/NewtonRingsI.m | 64 | ||||
| -rw-r--r-- | transverse/SimpleGaussian.m | 54 | ||||
| -rw-r--r-- | transverse/SimpleGaussian_rdep.m | 82 | ||||
| -rw-r--r-- | transverse/complexLaguerreGaussianE.m | 197 | ||||
| -rw-r--r-- | transverse/decompose.m | 174 | ||||
| -rw-r--r-- | transverse/fitgaussianbeamrad.m | 19 | ||||
| -rw-r--r-- | transverse/hermitepoly.m | 18 | ||||
| -rw-r--r-- | transverse/overlap.m | 83 | ||||
| -rw-r--r-- | transverse/polarmesh.m | 92 | ||||
| -rw-r--r-- | transverse/recompose.m | 83 |
16 files changed, 1345 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);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/transverse/AiryI.m b/transverse/AiryI.m new file mode 100644 index 0000000..d3eb610 --- /dev/null +++ b/transverse/AiryI.m @@ -0,0 +1,25 @@ +%---------------------------------------------------------------
+% PROGRAM: AiryI
+% AUTHOR: Andri M. Gretarsson
+% DATE: 8/7/03
+%
+%
+% SYNTAX: I=AiryI([a,w0,L,lambda <,N>],r);
+% <,...> indicates optional argument
+%
+% Returns the intensity of an Airy disk as a function of radius
+%
+% 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: I=AiryI([a,w0,L,lambda <,N>],r);
+%---------------------------------------------------------------
+
+function I=AiryI(params,r);
+
+I=AiryE(params,r).*conj(AiryE(params,r));
\ No newline at end of file diff --git a/transverse/AiryI_Fraunhofer.m b/transverse/AiryI_Fraunhofer.m new file mode 100644 index 0000000..d0a55e2 --- /dev/null +++ b/transverse/AiryI_Fraunhofer.m @@ -0,0 +1,33 @@ +%---------------------------------------------------------------
+% PROGRAM: AiryI_Fraunhofer
+% AUTHOR: Andri M. Gretarsson
+% DATE: 1/30/06
+%
+%
+% SYNTAX: I=AiryI([a,L,lambda],r);
+% <,...> indicates optional argument
+%
+% Returns the intensity of an Airy disk as a function of radius
+% in the Fraunhofer (small aperture) and paraxial (small angle)
+% limits.
+%
+% a = radius of aperture
+% L = propagation distance (axially) from aperture
+% lambda = wavelength
+%
+% Last updated: 1/30/06 by AMG
+%---------------------------------------------------------------
+%% SYNTAX: I=AiryI_Fraunhofer([a,L,lambda],r);
+%---------------------------------------------------------------
+
+function I=AiryI(params,r);
+
+a=params(1);
+L=params(2);
+lambda=params(3);
+k=2*pi/lambda;
+
+r=(r==0)*lambda/1000+r; %Gives the limit of Bessel(1,r)/r as r->0 (can't evaluate BesselJ(1,0)/0).
+I=(2*BesselJ(1,k*a*r/L)./(k*a*r./L)).^2;
+
+%I is normalized to the intensity in the center.
\ No newline at end of file diff --git a/transverse/HermiteGaussianE.m b/transverse/HermiteGaussianE.m new file mode 100644 index 0000000..ced8f6f --- /dev/null +++ b/transverse/HermiteGaussianE.m @@ -0,0 +1,134 @@ +function z=HermiteGaussianE(params,x,varargin)
+% Returns the complex field amplitude of a single polarization Hermite-Gaussian mode.
+%----------------------------------------------------------------------------------------------
+% PROGRAM: HermiteGaussianE
+% AUTHOR: Andri M. Gretarsson
+% DATE: 6/26/03
+%
+% SYNTAX: z=HermiteGaussianE([l,m,q <,lambda,a>],x <,y>);
+% <...> indicates optional arguments
+%
+% Returns the complex field amplitude of a single polarization Hermite-Gaussian mode.
+% The domain of the function is specified in cartesian coordinates x,y.
+% If the input argument a is not specified or is equal to 1, the area under the
+% surface z(x,y) is equal to 1. If z=z(x), i.e. the one dimensional form is being used,
+% then the area under the surface formed by rotating z(x) around the z axis equals 1.
+% In short the functions are normalized and form an orthonormal basis wrt. the inner
+% product formed by multiplying two HermiteE's together and integrating over area.
+% The formula is taken from "Principles of Lasers" by Orazio Svelto, 4th ed. Section 5.5.1.3.
+%
+% Note that it is possible to specify the parameters l,m,q,a,lambda as column vectors. In that case,
+% the matrix returned z, will be a stack of 2D matrixes (i.e. a 3D matrix), where successive 2D matrixes
+% in the stack (i.e. successive layers of the 3D matrix) correpsond to successive rows of the
+% l,m,q,a,lambda columns
+%
+% INPUT ARGUMENTS:
+% ----------------
+% l,m = Hermite Gaussian mode numbers (integers, positive or zero). Can be a scalar or a column vector.
+% q = complex radius of curvature of the beam (1/q = 1/R + i*lambda/pi/w^2).
+% Can be a scalar or a column vector.
+% a = complex amplitude (includes phase, e.g. for a beam that
+% has been propageted with an ABCD matrix a = 1/(A + [B/q1])^(1+l+m) ).
+% See for example, O. Svelto, Principles of Lasers 4th ed. eqn. 4.7.30.
+% Can be a scalar or a column vector.
+% lambda = Wavelength of the light. Can be a scalar or a column vector. Default is 1064e-9;
+% x = x position vector or a 2D matrix generated by meshgrid.
+% y = y position vector or a 2D matrix generated by meshgrid.
+% If y is not specified then the default y is used:
+% if x is 1D, y=zeros(size(x)). If x is 2D, y=x.
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% z(i,j) = Complex field of the Hermite Gaussian mode with parameters l,m,... at
+% ( x(i,j),y(i,j) ). May be a vector or a matrix depending on whether
+% x and y are vectors or matrixes. If in addition, l,m,q,... are columns,
+% then z(i,j,k) will correspond to the complex field amplitude of the
+% Hermite Gaussian mode corresponding to the parameters l(k),m(k),...
+%
+% NOTES:
+% ------
+% If x and y are not vectors but matrixes generated by
+% the matlab function meshgrid, then the output variable z
+% is a matrix rather than a vector. The matrix form allows
+% the function HermiteGaussianE to have a plane as it's domain
+% rather than a curve.
+%
+% If the parameters l,m,q,p,lambda are equal length column vectors rather
+% than scalars, z is a size(x)*length(lambda) matrix. E.g. if size(x) is n*n
+% then each level z(:,:,k) is a 2D field of a Hermite Gaussian with
+% the parameters given by [l(k),m(k),q(k),lambda(k),a(k)].
+%
+% EXAMPLE 1 (2D : two TEM modes):
+% w=[0.001; 0.001];
+% xseed=[-3*max(w):max(w)/30:3*max(w)]; yseed=xseed';
+% [x,y]=meshgrid(xseed,yseed);
+% lambda = [1.064e-6 ; 1.064e-6];
+% R = [-30 ; -30];
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1;1];
+% l=[1;2]; m=[0;2];
+% E=HermiteGaussianE([l,m,q,lambda,a],x,y);
+% subplot(2,1,1); h1=pcolor(x,y,abs(E(:,:,1))); set(h1,'EdgeColor','none'); axis square;
+% subplot(2,1,2); h2=pcolor(x,y,abs(E(:,:,2))); set(h2,'EdgeColor','none'); axis square; shg;
+%
+% EXAMPLE 2 (1D : four TEM modes):
+% w=[1,2,3,4].'; x=[-10:0.001:10]; lambda=[1,1,1,1].'*656e-9; R=[1,1,1,1].'*1000;
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1,1,1,1].'; l = [0,1,2,3].'; m=[0,0,0,0].';
+% E=HermiteGaussianE([l,m,q,lambda,a],x); I=E.*conj(E); phi=angle(E);
+% plot(x,I(:,1),x,I(:,2),x,I(:,3),x,I(:,4));
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+%% SYNTAX: z=HermiteGaussianE([l,m,q <,lambda,a>],x <,y>);
+%----------------------------------------------------------------------------------------------
+
+defaultcoord2=0;
+if nargin>=3, y=varargin{1}; else defaultcoord2=1; end
+
+if min(size(x))==1
+ if size(x,1)<size(x,2), x=x'; end %make x and y columnar
+ if defaultcoord2
+ y=zeros(size(x));
+ else
+ if size(y,1)<size(y,2), y=y'; end
+ end
+
+end
+
+if min(size(x)) > 1
+ z=zeros(size(x,1),size(x,2),size(params,1)); % need this since zeros(size(y),10) gives a 2D matrix even if y is 2D! (Matlab feature.)
+ if defaultcoord2, y=transpose(x); end
+else
+ z=zeros(length(x),size(params,1));
+end
+
+l=params(:,1);
+m=params(:,2);
+q=params(:,3);
+if size(params,2)>=4
+ lambda=params(:,4);
+else
+ lambda=1064e-9;
+end
+if size(params,2)>=5
+ a=params(:,5);
+else
+ a=ones(size(q));
+end
+
+w=w_(q,lambda);
+
+if min(size(x))>=2 % x is a matrix (2D array)
+ for u=1:size(params,1)
+ z(:,:,u) = a(u)...
+ .* sqrt((2/pi))/2^((l(u)+m(u))/2)/sqrt(factorial(l(u))*factorial(m(u)))/w(u)...
+ .* hermitepoly(l(u),sqrt(2).*x/w(u)).*hermitepoly(m(u),sqrt(2).*y/w(u))...
+ .* exp( -i*2*pi/lambda(u)*(x.^2+y.^2)/2/q(u) ) ;
+ end
+else
+ for u=1:size(params,1) % x is a vector (1D array)
+ z(:,u) = a(u)...
+ .* sqrt((2/pi))/2^((l(u)+m(u))/2)/sqrt(factorial(l(u))*factorial(m(u)))/w(u)...
+ .* hermitepoly(l(u),sqrt(2).*x/w(u)).*hermitepoly(m(u),sqrt(2).*y/w(u))...
+ .* exp( -i*2*pi/lambda(u)*(x.^2+y.^2)/2/q(u) );
+ end
+end
diff --git a/transverse/LaguerreGaussianE.m b/transverse/LaguerreGaussianE.m new file mode 100644 index 0000000..f95985d --- /dev/null +++ b/transverse/LaguerreGaussianE.m @@ -0,0 +1,189 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: LaguerreGaussianE
+% AUTHOR: Andri M. Gretarsson
+% DATE: 6/26/03
+%
+% SYNTAX: z=LaguerreGaussianE([p,m,q <,lambda,a>],r <,theta,coordtype>);
+% <...> indicates optional arguments
+%
+% Returns the complex field amplitude of a Laguerre-Gaussian mode as a function
+% of polar coordinates r and theta. Formula adapted from A. E. Siegman,
+% "Lasers" 1st ed. eqn. 64 in Ch. 16. I leave out the gouy phase
+% factor since it is meaningless except as a relative phase difference
+% between axially separated parts of the same beam. In other words it
+% only appears when propagating the beam. The function prop.m does
+% both the q transformation and supplies the appropriate gouy phase.
+% This factor and other phase and amplitude factors can be included via the
+% complex argument a if desired. Finally, this function can also be called
+% with cartesian coordinates.
+%
+%
+% INPUT ARGUMENTS:
+% ----------------
+% p,m = Laguerre Gaussian mode numbers (column vector)
+% q = complex radius of curvature of the beam (column vector)
+% lambda= Wavelength of the light (column vector)
+% a = complex prefactor ( includes gouy phase, e.g. for a LG beam that
+% has been propageted with an ABCD matrix a = 1/(A + [B/q1])^(1+2p+m) ).
+% (column vector).
+% r = radial position vector (or matrix generated by meshgrid)
+% theta = azimuthal position vector (or matrix generated by meshgrid).
+% If r is 1D and theta is not specified the default theta is used,
+% theta=zeros(size(r)). If r is a 2D mesh theta must also be specified.
+% The function polarmesh.m can be usefule for generating the r and
+% theta meshes.
+% coordtype = (string) label for the type of coordinates supplied in r and theta. If
+% this argument is not specified or is equal to 'pol', r and theta are assumed
+% to be polar coordinates (r,theta) as described above. If on the other
+% hand coordtype='cart' then r is assumed to be the cartesian x coordinate
+% and theta the cartesian y coordinate. This can be useful e.g. if one wants to
+% specify an x or y shift in the position of the mode. (See for example
+% Laguerre_demo.m). coordtype is normally the 8th input argument but can
+% be specified as the 7th input argument instead if the default y is desired
+% (y=transpose(x)).
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% z(i,j)= Complex field of the Laguerre Gaussian mode at
+% ( r(i,j),theta(i,j) ). May be a vector or
+% a matrix depending on whether r and theta
+% are vectors or matrixes.
+%
+%
+% NOTES:
+% ------
+% If r and theta are not vectors but matrixes generated by
+% the matlab function meshgrid, then the output variable z
+% is a matrix rather than a vector. The matrix form allows
+% the function HermiteGaussianE to have a plane as it's domain
+% rather than a curve.
+%
+% If the parameters p,m,q,p,lambda are equal length column vectors rather
+% than scalars z is a size(r)*length(lambda) matrix. E.g. if size(r) is n*n
+% then each level z(:,:,k) is a 2D field of a Laguerre Gaussian with
+% the parameters given by [l(k),m(k),q(k),lambda(k),a(k)].
+%
+% EXAMPLE 1 (2D, polar domain):
+% w=[0.001; 0.001];
+% rseed=[0*max(w):max(w)/30:3*max(w)]; thetaseed=[0:360]*pi/180;
+% [r,theta]=meshgrid(rseed,thetaseed);
+% lambda = [1.064e-6 ; 1.064e-6];
+% R = [-30 ; -30];
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1;1];
+% p=[1;2]; m=[0;2];
+% E=LaguerreGaussianE([p,m,q,lambda,a],r,theta)+LaguerreGaussianE([p,-m,q,lambda,a],r,theta);
+% [x,y]=pol2cart(theta,r); colormap(bone);
+% subplot(2,1,1); h1=pcolor(x,y,abs(E(:,:,1))); set(h1,'EdgeColor','none'); axis square;
+% subplot(2,1,2); h2=pcolor(x,y,abs(E(:,:,2))); set(h2,'EdgeColor','none'); axis square; shg;
+%
+% EXAMPLE 2 (1D, cartesian domain, default y):
+% w=[1,2,3,4].'; x=[-10:0.001:10]; lambda=[1,1,1,1].'*656e-9; R=[1,1,1,1].'*1000;
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1,1,1,1].'; p = [0,1,2,3].'; m=[0,0,0,0].';
+% E=LaguerreGaussianE([p,m,q,lambda,a],x,'cart'); I=E.*conj(E); phi=angle(E);
+% plot(x,I(:,1),x,I(:,2),x,I(:,3),x,I(:,4));
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+%% SYNTAX: z=LaguerreGaussianE([p,m,q <,lambda,a>],r <,theta,coordtype>);
+%----------------------------------------------------------------------------------------------
+
+function z=LaguerreGaussianE(params,r,varargin);
+
+if nargin>=3
+ if isstr(varargin{1})
+ if strcmp(varargin{1},'cart') % use cartesian coordinates
+ defaultcoord2=1;
+ cartesianflag=1;
+ else % use polar coordinates with the default theta
+ defaultcoord2=1;
+ cartesianflag=0;
+ end
+ else % use polar coordinates with the specified theta
+ defaultcoord2=0;
+ cartesianflag=0;
+ theta=varargin{1};
+ end
+else % use polar coordinates with the default theta
+ defaultcoord2=1;
+ cartesianflag=0;
+end
+
+if nargin>=4
+ defautlcoord2=0;
+ if isstr(varargin{2}) & strcmp(varargin{2},'cart')
+ cartesianflag=1; % use cartesian coordinates with specified y
+ else % use polar coordinates with the specified theta
+ cartesianflag=0;
+ end
+end
+
+
+if cartesianflag % cartesian (x,y) domain supplied
+ x=r;
+ if min(size(x))==1 % map is 2->1 on a cartesian domain
+ if size(x,1)<size(x,2), x=x'; end % make x and y columnar
+ if defaultcoord2
+ y=zeros(size(x));
+ else
+ y=theta;
+ if size(y,1)<size(y,2), y=y'; end
+ end
+ end
+ if min(size(x)) > 1 % map is 2->2 on a cartesian domain
+ if defaultcoord2
+ y=transpose(x);
+ else
+ y=theta;
+ end
+ z=zeros(size(x,1),size(x,2),size(params,1)); % need this since zeros(size(y),10) gives a 2D matrix even if y is 2D! (Matlab feature.)
+ else
+ z=zeros(size(x),size(params,1));
+ end
+ [theta,r]=cart2pol(x,y); % convert to polar coords for calculation
+else % polar (r,theta) domain supplied
+ if min(size(r))==1 % map is 2->1 on a polar domain
+ if size(r,1)<size(r,2), r=r.'; end % make r columnar
+ if defaultcoord2 % make theta columnar
+ theta=zeros(size(r)); % default 1D theta is zero
+ else
+ if size(theta,1)<size(theta,2), theta=theta.'; end
+ end
+ else % otherwise assume r and theta are already in meshgrid format
+ z=zeros(size(r,1),size(r,2),size(params,1)); % need this since zeros(size(r),10) gives a 2D matrix even if y is 2D! (Matlab feature.)
+ end
+end
+
+
+p=params(:,1);
+m=params(:,2);
+signm=sign(m);
+m=abs(m);
+q=params(:,3);
+if size(params,2)>=4
+ lambda=params(:,4);
+else
+ lambda=1064e-9;
+end
+if size(params,2)>=5
+ a=params(:,5);
+else
+ a=ones(size(q));
+end
+
+w=w_(q,lambda);
+
+if min(size(r))>=2
+ for u=1:size(params,1)
+ z(:,:,u) = a(u)...
+ .* sqrt(2*factorial(p(u))/(1+(m(u)==0))/pi/(factorial( m(u)+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^m(u) .*exp(1i*signm(u)*m(u).*theta).* LaguerrePoly([p(u),m(u)],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u));
+ end
+else
+ for u=1:size(params,1)
+ z(:,u) = a(u)...
+ .* sqrt(2*factorial(p(u))/(1+(m(u)==0))/pi/(factorial( m(u)+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^m(u) .* exp(i*signm(u)*m(u).*theta).* LaguerrePoly([p(u),m(u)],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u));
+ end
+end
diff --git a/transverse/LaguerrePoly.m b/transverse/LaguerrePoly.m new file mode 100644 index 0000000..b29d3d0 --- /dev/null +++ b/transverse/LaguerrePoly.m @@ -0,0 +1,24 @@ +%-----------------------------------
+% Associated Laguerre Polynomial
+% (See e.g. Arfken section 13.2)
+%
+% SYNTAX y=Lnk([n,k],x)
+%-----------------------------------
+
+function y=LaguerrePoly(params,x)
+
+n=params(1);
+k=params(2);
+
+m=[0:n];
+
+a=factorial(n+k)*ones(1,length(m));
+b=factorial(n-m);
+c=factorial(k+m);
+d=factorial(m);
+e=(-1).^m;
+
+y=zeros(size(x));
+for s=1:n+1;
+ y = y + a(s) ./ b(s) ./ c(s) ./ d(s) .* e(s) .* x.^m(s);
+end
diff --git a/transverse/NewtonRingsI.m b/transverse/NewtonRingsI.m new file mode 100644 index 0000000..961495b --- /dev/null +++ b/transverse/NewtonRingsI.m @@ -0,0 +1,64 @@ +%---------------------------------------------------------------
+% PROGRAM: NRI
+% AUTHOR: Andri M. Gretarsson
+% DATE: 8/5/03
+%
+%
+% SYNTAX: z=NwetonRingsI([w1,w2,R1,R2,P1,P2,lambda,deltaphi],r);
+%
+% Returns the intensity of a set of Newton's rings as a function
+% of the radial coordinate R. (See e.g. Fundamentals of Physical Optics
+% by Jenkins and White, 1st ed. Section 4.4). The Newton's rings are
+% generated by the interference of two coaxial beams with different radii
+% of curvature and/or width.
+%
+% w1 = 1/e field amplitude radius of the 1st of the interfering beams
+% w2 = 1/e field amplitude radius of the 2nd of the interfering beams
+% R1 = Radius of curvature of the phasefront of the 1st of the interfering beams
+% R2 = Radius of curvature of the phasefront of the 2nd of the interfering beams
+% P1 = Radius of curvature of the phasefront of the 2nd of the interfering beams
+% P2 = Radius of curvature of the phasefront of the 2nd of the interfering beams
+% lambda = Radius of curvature of the phasefront of the 2nd of the interfering beams
+% deltaphi = Radius of curvature of the phasefront of the 2nd of the interfering beams
+% z(k) = Intensity of the interference pattern at radiusr(k).
+%
+% If x and y are not vectors but matrixes generated by
+% the matlab function meshgrid, then the output variable z
+% is a matrix rather than a vector. The matrix form allows
+% the function to have a plane as it's domain rather than
+% a curve. To generate a plane domain one uses meshgrid.
+%
+% EXAMPLE: thetaoffset=20*pi/180;
+% rseed=[0:0.01:sqrt(0.3)].^2;
+% thetaseed=[0:1:360]*pi/180;
+% [theta,r]=meshgrid(thetaseed,rseed);
+% [x,y]=pol2cart(theta,r);
+% h=pcolor(x,y,NewtonRingsI([0.01,100,2e3,1e99,1,1,1.024e-6,thetaoffset],r));
+% colormap('bone');
+% set(h,'EdgeColor','none');
+% set(h,'FaceColor','interp');
+% set(gca,'Visible','off');
+% set(gcf,'Color','black');
+% axis square
+% shg;
+%
+% Last updated: 8/5/03 by AMG
+%
+%---------------------------------------------------------------
+%% SYNTAX: z=NewtonRingsI([w1,w2,R1,R2,P1,P2,lambda,deltaphi],r);
+%---------------------------------------------------------------
+
+function z=NewtonRingsI(params,r)
+
+w1=params(1);
+w2=params(2);
+R1=params(3);
+R2=params(4);
+P1=params(5);
+P2=params(6);
+lambda=params(7);
+deltaphi=params(8);
+a=1/2/R1-1/2/R2;
+
+z = exp(-r.^2/w1).*exp(-r.^2/w2).*...
+ ( sin( (a/2/lambda)* r.^2 - deltaphi) ).^2 + abs(P1-P2);
\ No newline at end of file diff --git a/transverse/SimpleGaussian.m b/transverse/SimpleGaussian.m new file mode 100644 index 0000000..327f11e --- /dev/null +++ b/transverse/SimpleGaussian.m @@ -0,0 +1,54 @@ +%---------------------------------------------------------------
+% Field of the TEM00 Gaussian mode as a function of distance
+% from waist and axial distance. (See e.g. Orazio Svelto,
+% Principles of Lasers, 4th ed. page 152, Eqn's 4.715a-4.717c)
+%
+% SYNTAX: [E,w,R,phi,zR]=SimpleGaussian([w0,lambda],z,r);
+%
+% INPUT ARGUMENTS:
+% w0 = Gaussian field radius at waist
+% z = axial distance from waist
+% lambda = wavelength
+%
+% z = distance from waist (Nx1 vector)
+% r = distance from beam axis (Nx1 vector)
+%
+% OUTPUT ARGUMENTS:
+% E = complex electric field normalized to the field
+% amplitude at the center of the waist
+% w = width of the beam (radius at which the field amplitude
+% falls to 1/e of it's value on the beam axis
+% R = Radius of curvature of phasefront
+% phi = Guoy phase
+% zR = Raleigh range
+%
+%---------------------------------------------------------------
+% SYNTAX: [E,w,R,phi,zR]=SimpleGaussian([w0,lambda],z,r);
+%---------------------------------------------------------------
+
+function [E,w,R,phi,zR]=SimpleGaussian(params,z,r);
+
+w0=params(1); %Beam field width at waist
+lambda=params(2); %Wavelength
+
+higherorder=0;
+if length(params)>=3,
+ l=params(3);
+ m=params(4);
+ higherorder=1;
+end
+
+E=zeros(length(z),length(r));
+w=zeros(length(z),1);
+R=zeros(length(z),1);
+phi=zeros(length(z),1);
+
+if higherorder==0
+ for s=1:length(z)
+ [E(s,:),w(s),R(s),phi(s),zR]=SimpleGaussian_rdep([w0,z(s),lambda],r);
+ end
+else
+ for s=1:length(z)
+ [E(s,:),w(s),R(s),phi(s),zR]=SimpleGaussian_rdep([w0,z(s),lambda,0,l,m],r);
+ end
+end
\ No newline at end of file diff --git a/transverse/SimpleGaussian_rdep.m b/transverse/SimpleGaussian_rdep.m new file mode 100644 index 0000000..ab81f21 --- /dev/null +++ b/transverse/SimpleGaussian_rdep.m @@ -0,0 +1,82 @@ +%---------------------------------------------------------------
+% Field of the TEMnm Gaussian mode as a function of radial
+% distance from beam axis.
+%
+% SYNTAX: E=SimpleGaussian_rdep([w0,z,lambda,n,m],x,y);
+%
+% INPUT ARGUMENTS:
+% w0 = Gaussian field radius at waist
+% z = distance from beam axis
+% lambda = wavelength
+%
+% r = distance from beam axis (Nx1 vector)
+% n = horizontal mode number
+% m = vertical mode number
+%
+% OUTPUT ARGUMENTS:
+% E = complex electric field normalized to the field
+% amplitude at the center of the waist
+% w = width of the beam (radius at which the field amplitude
+% falls to 1/e of it's value on the beam axis
+% R = Radius of curvature of phasefront
+% phi = Guoy phase
+% zR = Raleigh range
+%
+%---------------------------------------------------------------
+% SYNTAX: E=SimpleGaussian_rdep([w0,z,lambda,twoD,n,m],x,y);
+%---------------------------------------------------------------
+
+function [E,w,R,phi,zR]=SimpleGaussian_rdep(params,x,varargin);
+
+if abs(params(2))<=1e-99, params(2)=1e-99; end
+ %Prevents divide by zero error at waist
+
+if nargin>=3
+ y=varargin{1};
+else
+ y=zeros(size(x));
+end
+higherorder=0;
+
+w0=params(1); %Beam field width at waist
+z=params(2); %Distance from beam axis
+if length(params)>=4
+ twoD=params(4);
+else
+ twoD=0;
+end
+if length(params)>=5 %Will be calculating higher order mode
+ l=params(5);
+ m=params(6);
+ higherorder=1;
+end
+
+
+lambda=params(3); %Wavelength
+k=2*pi/lambda; %Wavenumber
+zR=pi*w0^2/lambda; %Raleigh range
+
+w=w0*sqrt(1+(z/zR)^2); %Beam (field) width
+R=z*(1+(zR/z)^2); %Radius of curvature
+phi=atan(z/zR); %Guoy phase
+
+if twoD~=1
+ if higherorder~=1
+ E = (w0/w) * exp(-(x.^2+y.^2)/w^2) .* exp(-i*k*(x.^2+y.^2)/2/R) * exp( i*phi);
+ else
+ E = (w0/w) * hermitepoly(l,sqrt(2)*x/w) .* hermitepoly(m,sqrt(2)*y/w)...
+ .*exp(-(x.^2+y.^2)/w^2) .* exp(-i*k*(x.^2+y.^2)/2/R) * exp(i*(1+l+m)*phi);
+ end
+else
+ E=zeros(length(y),length(x));
+ if higherorder~=1
+ for s=1:length(y)
+ E(s,:) = (w0/w) * exp(-(x.^2+y(s).^2)/w^2) .* exp(-i*k*(x.^2+y(s).^2)/2/R) * exp( i*phi);
+ end
+ else
+ for s=1:length(y)
+ E(s,:) = (w0/w) * hermitepoly(l,sqrt(2)*x/w) .* hermitepoly(m,sqrt(2)*y(s)/w)...
+ .*exp(-(x.^2+y(s).^2)/w^2) .* exp(-i*k*(x.^2+y(s).^2)/2/R) * exp(i*(1+l+m)*phi);
+ end
+ end
+end
diff --git a/transverse/complexLaguerreGaussianE.m b/transverse/complexLaguerreGaussianE.m new file mode 100644 index 0000000..e6dd73e --- /dev/null +++ b/transverse/complexLaguerreGaussianE.m @@ -0,0 +1,197 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: LaguerreGaussianE
+% AUTHOR: Andri M. Gretarsson
+% DATE: 6/26/03
+%
+% SYNTAX: z=LaguerreGaussianE([p,m,q <,lambda,a>],r <,theta,coordtype>);
+% <...> indicates optional arguments
+%
+% Returns the complex field amplitude of a Laguerre-Gaussian mode as a function
+% of polar coordinates r and theta. Formula adapted from A. E. Siegman,
+% "Lasers" 1st ed. eqn. 64 in Ch. 16. I leave out the gouy phase
+% factor since it is meaningless except as a relative phase difference
+% between axially separated parts of the same beam. In other words it
+% only appears when propagating the beam. The function prop.m does
+% both the q transformation and supplies the appropriate gouy phase.
+% This factor and other phase and amplitude factors can be included via the
+% complex argument a if desired. Finally, this function can also be called
+% with cartesian coordinates.
+%
+%
+% INPUT ARGUMENTS:
+% ----------------
+% p,m = Laguerre Gaussian mode numbers (column vector)
+% q = complex radius of curvature of the beam (column vector)
+% lambda= Wavelength of the light (column vector)
+% a = complex prefactor ( includes gouy phase, e.g. for a LG beam that
+% has been propageted with an ABCD matrix a = 1/(A + [B/q1])^(1+2p+m) ).
+% (column vector).
+% r = radial position vector (or matrix generated by meshgrid)
+% theta = azimuthal position vector (or matrix generated by meshgrid).
+% If r is 1D and theta is not specified the default theta is used,
+% theta=zeros(size(r)). If r is a 2D mesh theta must also be specified.
+% The function polarmesh.m can be usefule for generating the r and
+% theta meshes.
+% coordtype = (string) label for the type of coordinates supplied in r and theta. If
+% this argument is not specified or is equal to 'pol', r and theta are assumed
+% to be polar coordinates (r,theta) as described above. If on the other
+% hand coordtype='cart' then r is assumed to be the cartesian x coordinate
+% and theta the cartesian y coordinate. This can be useful e.g. if one wants to
+% specify an x or y shift in the position of the mode. (See for example
+% Laguerre_demo.m). coordtype is normally the 8th input argument but can
+% be specified as the 7th input argument instead if the default y is desired
+% (y=transpose(x)).
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% z(i,j)= Complex field of the Laguerre Gaussian mode at
+% ( r(i,j),theta(i,j) ). May be a vector or
+% a matrix depending on whether r and theta
+% are vectors or matrixes.
+%
+%
+% NOTES:
+% ------
+% If r and theta are not vectors but matrixes generated by
+% the matlab function meshgrid, then the output variable z
+% is a matrix rather than a vector. The matrix form allows
+% the function HermiteGaussianE to have a plane as it's domain
+% rather than a curve.
+%
+% If the parameters p,m,q,p,lambda are equal length column vectors rather
+% than scalars z is a size(r)*length(lambda) matrix. E.g. if size(r) is n*n
+% then each level z(:,:,k) is a 2D field of a Laguerre Gaussian with
+% the parameters given by [l(k),m(k),q(k),lambda(k),a(k)].
+%
+% EXAMPLE 1 (2D, polar domain):
+% w=[0.001; 0.001];
+% rseed=[0*max(w):max(w)/30:3*max(w)]; thetaseed=[0:360]*pi/180;
+% [r,theta]=meshgrid(rseed,thetaseed);
+% lambda = [1.064e-6 ; 1.064e-6];
+% R = [-30 ; -30];
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1;1];
+% p=[1;2]; m=[0;2];
+% E=LaguerreGaussianE([p,m,q,lambda,a],r,theta)+LaguerreGaussianE([p,-m,q,lambda,a],r,theta);
+% [x,y]=pol2cart(theta,r); colormap(bone);
+% subplot(2,1,1); h1=pcolor(x,y,abs(E(:,:,1))); set(h1,'EdgeColor','none'); axis square;
+% subplot(2,1,2); h2=pcolor(x,y,abs(E(:,:,2))); set(h2,'EdgeColor','none'); axis square; shg;
+%
+% EXAMPLE 2 (1D, cartesian domain, default y):
+% w=[1,2,3,4].'; x=[-10:0.001:10]; lambda=[1,1,1,1].'*656e-9; R=[1,1,1,1].'*1000;
+% q = (1./R - i* lambda./pi./w.^2).^(-1); a=[1,1,1,1].'; p = [0,1,2,3].'; m=[0,0,0,0].';
+% E=LaguerreGaussianE([p,m,q,lambda,a],x,'cart'); I=E.*conj(E); phi=angle(E);
+% plot(x,I(:,1),x,I(:,2),x,I(:,3),x,I(:,4));
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+%% SYNTAX: z=LaguerreGaussianE([p,m,q <,lambda,a>],r <,theta,coordtype>);
+%----------------------------------------------------------------------------------------------
+
+function z=LaguerreGaussianE(params,r,varargin);
+
+if nargin>=3
+ if isstr(varargin{1})
+ if strcmp(varargin{1},'cart') % use cartesian coordinates
+ defaultcoord2=1;
+ cartesianflag=1;
+ else % use polar coordinates with the default theta
+ defaultcoord2=1;
+ cartesianflag=0;
+ end
+ else % use polar coordinates with the specified theta
+ defaultcoord2=0;
+ cartesianflag=0;
+ theta=varargin{1};
+ end
+else % use polar coordinates with the default theta
+ defaultcoord2=1;
+ cartesianflag=0;
+end
+
+if nargin>=4
+ defautlcoord2=0;
+ if isstr(varargin{2}) & strcmp(varargin{2},'cart')
+ cartesianflag=1; % use cartesian coordinates with specified y
+ else % use polar coordinates with the specified theta
+ cartesianflag=0;
+ end
+end
+
+
+if cartesianflag % cartesian (x,y) domain supplied
+ x=r;
+ if min(size(x))==1 % map is 2->1 on a cartesian domain
+ if size(x,1)<size(x,2), x=x'; end % make x and y columnar
+ if defaultcoord2
+ y=zeros(size(x));
+ else
+ y=theta;
+ if size(y,1)<size(y,2), y=y'; end
+ end
+ end
+ if min(size(x)) > 1 % map is 2->2 on a cartesian domain
+ if defaultcoord2
+ y=transpose(x);
+ else
+ y=theta;
+ end
+ z=zeros(size(x,1),size(x,2),size(params,1)); % need this since zeros(size(y),10) gives a 2D matrix even if y is 2D! (Matlab feature.)
+ else
+ z=zeros(size(x),size(params,1));
+ end
+ [theta,r]=cart2pol(x,y); % convert to polar coords for calculation
+else % polar (r,theta) domain supplied
+ if min(size(r))==1 % map is 2->1 on a polar domain
+ if size(r,1)<size(r,2), r=r.'; end % make r columnar
+ if defaultcoord2 % make theta columnar
+ theta=zeros(size(r)); % default 1D theta is zero
+ else
+ if size(theta,1)<size(theta,2), theta=theta.'; end
+ end
+ else % otherwise assume r and theta are already in meshgrid format
+ z=zeros(size(r,1),size(r,2),size(params,1)); % need this since zeros(size(r),10) gives a 2D matrix even if y is 2D! (Matlab feature.)
+ end
+end
+
+
+p=params(:,1);
+m=params(:,2);
+signm=sign(m);
+m=abs(m);
+q=params(:,3);
+if size(params,2)>=4
+ lambda=params(:,4);
+else
+ lambda=1064e-9;
+end
+if size(params,2)>=5
+ a=params(:,5);
+else
+ a=ones(size(q));
+end
+
+w=w_(q,lambda);
+
+if min(size(r))>=2
+ for u=1:size(params,1)
+ z(:,:,u) = a(u)...
+ .* sqrt(2*factorial(p(u))/(1+(m(u)==0))/pi/(factorial( m(u)+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^m(u) .*exp(1i*signm(u)*m(u).*theta).* LaguerrePoly([p(u),m(u)],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u))
+ + a(u)...
+ .* sqrt(2*factorial(p(u))/(1+((-m(u))==0))/pi/(factorial( (-m(u))+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^(-m(u)) .*exp(1i*signm(u)*(-m(u)).*theta).* LaguerrePoly([p(u),(-m(u))],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u));
+ end
+else
+ for u=1:size(params,1)
+ z(:,u) = a(u)...
+ .* sqrt(2*factorial(p(u))/(1+(m(u)==0))/pi/(factorial( m(u)+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^m(u) .*exp(1i*signm(u)*m(u).*theta).* LaguerrePoly([p(u),m(u)],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u))
+ + a(u)...
+ .* sqrt(2*factorial(p(u))/(1+((-m(u))==0))/pi/(factorial( (-m(u))+p(u) )))/w(u)...
+ .* (sqrt(2)*r/w(u)).^(-m(u)) .*exp(1i*signm(u)*(-m(u)).*theta).* LaguerrePoly([p(u),(-m(u))],2*r.^2/w(u).^2)...
+ .* exp( -1i*2*pi/lambda(u)*r.^2/2/q(u));
+ end
+end
diff --git a/transverse/decompose.m b/transverse/decompose.m new file mode 100644 index 0000000..e9b583a --- /dev/null +++ b/transverse/decompose.m @@ -0,0 +1,174 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: decompose
+% AUTHOR: Andri M. Gretarsson
+% DATE: 7/10/04
+%
+% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>]);
+% <...> indicates optional arguments
+%
+% Calculates the coefficients of the decomposition of the function z1 into
+% Hermite Gaussian modes or Laguerre Gaussian modes.
+%
+% INPUT ARGUMENTS:
+% ----------------
+% z1 = The values of the function to be decomposed. 2D (nxn) matrix
+% domain = the domain values at which the values in z1 are specified.
+% The domain is a nxmx2 array where domain(:,:,1) and
+% domain(:,:,2) are the x and y meshes corresponding to the
+% values in z1. (These meshes are often generated using meshgrid.m.)
+% type = 'hg' for a Hermite Gaussian mode decomposition
+% 'lg' for a Laguerre Gaussian mode decomposition
+% terms = This argument can be specified either as a scalar or as a matrix.
+% If it is specified as a scalar, it indicates the number of
+% terms out to which to calculate the decomposition. Since the
+% Hermite and Laguerre bases are two dimensional and are specified
+% by a 2D argument (l,m) or (p,m), default rules are applied to
+% decide exactly which terms to calculate. The output
+% matrix tmat indicates what terms were calculated.
+%
+% When terms is specified as a matrix, it indicates exactly which
+% terms to calculate. For Hermite Gaussian modes (type='hg')
+% terms is a LxM matrix where a 1 in the (i,j)th entry indicates
+% that the term TEM(l,m)=(i-1,j-1) will be calculated. A zero
+% indicates it will not be calculated. For example:
+%
+% 1 0 0
+% terms = 0 1 0
+% 0 1 0
+%
+% indicates that the coefficients of the Hermite Gaussian modes
+% (0,0), (1,1) and (2,1) will be calculated but no others.
+%
+% For Laguerre Gaussian modes (type='lg') the terms matrix is
+% interpreted similarly. However, since the (p,m) mode can
+% have either positive or negative m terms the terms matrix is now
+% a PxMx2 matrix, where terms(:,:,1) indicates which positive-m
+% terms to calculate and terms(:,:,2) indicates which negative-m
+% terms to calculate. However, note that since m=0 must not
+% be duplicated, the first column in terms(:,:,2) is always zero.
+% Thus for example:
+%
+% 1 0 0
+% terms(:,:,1) = 0 1 0
+% 0 1 0
+%
+% 0 0 0
+% terms(:,:,2) = 0 1 0
+% 0 1 1
+%
+% indicates that the coefficients of the Laguerre Gaussian modes
+% (0,0), (1,1), (2,1), (1,-1), (2,-1) and (2,-2) will be calculated.
+%
+% q = the complex radius of curvature "q" of the Gaussian basis.
+% lambda = wavelength of the light in the Gaussian basis. Default is 1.064 microns
+% accuracy = only calculate the coefficients of the decomposition to the
+% specified accuracy. Actually rounds each result to the nearest
+% increment of accuracy. For example, if accuracy=0.3, then a coefficient
+% of 1.54 would be rounded to 1.5 while a coefficient of 1.56 would be
+% rounded to 1.8. Accuracy of 0 applies no rounding. Specifying
+% accuracy other than "0" does not speed up the calculation.
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% coeffs = the coefficients of the terms in the expansion. coeffs is a matrix
+% the same size a tmat.
+% tmat = coefficient request matrix. Often this would be specified by the user
+% by the input argument "terms", in other words tmat=terms. However, as explained
+% under the definition of the input argument terms, this function
+% allows the user to specify only how many terms he wants, in other words terms can be
+% a scalar and not a matrix. In this case tmat must be calculated according to
+% default rules. In this case it is convenient to have the terms request matrix
+% as an output argument. Each entry in the coefficient request matrix tmat correspond
+% to a particular term of the expansion (as explained above). If the tmat(i,j) is 1,
+% the coefficient of that term is calculated and returned in coeffs(i,j). If tmat(i,j)
+% is 0, the corresponding coefficient was not calculated and coeffs(i,j) is set to 0
+% regardless of whether the contribution of that term to z1 is really zero.
+%
+% EXAMPLE 1 (Hermite Gaussian, only number of terms specified):
+% [x,y]=meshgrid([-pi/2:0.1:pi/2],[-pi/2:0.1:pi/2]); z1=cos(sqrt(x.^2+y.^2));
+% clear domain; domain(:,:,1)=x; domain(:,:,2)=y;
+% w=pi/4; R=1e3; lambda=1e-6; q=(1./R - i* lambda./pi./w.^2).^(-1);
+% [coeffs,tmat]=decompose(z1,domain,'hg',140,[q,lambda,1e-6])
+% z1recomposed=recompose(domain,'hg',coeffs,[q,lambda,1e-6]);
+% subplot(3,1,1); h=pcolor(x,y,abs(z1).^2); set(h,'EdgeColor','none'); axis square; colorbar
+% subplot(3,1,2); h=pcolor(x,y,abs(z1recomposed).^2); set(h,'EdgeColor','none'); axis square; colorbar
+% subplot(3,1,3); h=pcolor(x,y,abs(z1).^2-abs(z1recomposed).^2); set(h,'EdgeColor','none'); axis square; colorbar; shg;
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>]);
+%----------------------------------------------------------------------------------------------
+function [coeffs,tmat]=decompose(z1,domain,type,terms,params)
+
+% HERMITE GAUSSIAN EXPANSION --------------------------------------------------------------------------
+if strcmpi(type,'hg')
+ q=params(1);
+ if length(params)>=2, lambda=params(2); else lambda=1.064e-6; end
+ if length(params)>=3, accuracy=params(3); else accuracy=1e-4; end
+
+ if size(terms,1)==1 && size(terms,2)==1 % Make terms request matrix
+ n=ceil(sqrt(terms));
+ tmat=ones(n,n);
+ if n^2>terms
+ tmat(end,end-ceil((n^2-terms)/2)+1:end)=0;
+ tmat(end-floor((n^2-terms)/2):end-1,end)=0;
+ end
+ else
+ tmat=terms;
+ end
+ coeffs=zeros(size(tmat));
+ for s=1:size(tmat,1) % Calculate the coefficients
+ l=s-1;
+ for t=1:size(tmat,2)
+ m=t-1;
+ if tmat(s,t)==1
+ z2=HermiteGaussianE([l,m,q,lambda],domain(:,:,1),domain(:,:,2));
+ coeffs(s,t)=overlap(z1,conj(z2),domain,1,accuracy);
+ else
+ coeffs(s,t)=0;
+ end
+ end
+ end
+
+% LAGUERRE GAUSSIAN EXPANSION --------------------------------------------------------------------------
+elseif strcmpi(type,'lg')
+ q=params(1);
+ if length(params)>=2, lambda=params(2); else lambda=1.064e-6; end
+ if length(params)>=3, accuracy=params(3); else accuracy=1e-4; end
+ if size(terms,1)==1 && size(terms,2)==1 % Make terms request matrix
+ n=ceil( (1+sqrt(1+8*terms))/4 );
+ tmat=ones(n,n,2); tmat(:,1,2)=0;
+ ndiff=2*n^2-n-terms;
+ if ndiff>0
+ rdiff=ceil(ndiff/2); ldiff=floor(ndiff/2);
+ tmat(end,end-ceil(rdiff/2)+1:end,1)=0;
+ tmat(end-floor(rdiff/2):end-1,end,1)=0;
+ tmat(end,end-ceil(ldiff/2)+1:end,2)=0;
+ tmat(end-floor(ldiff/2):end-1,end,2)=0;
+ end
+ %disp(' '); dispmat(tmat(:,:,1)); disp(' '); dispmat(tmat(:,:,2)); disp(' '); disp(num2str(sum(sum(sum(tmat)))));
+ else
+ tmat=terms;
+ end
+ clear coeffs;
+ coeffs(:,:,1)=zeros(size(tmat,1),size(tmat,2));
+ coeffs(:,:,2)=zeros(size(tmat,1),size(tmat,2));
+ for s=1:size(tmat,1) % Calculate the coefficients
+ p=s-1;
+ for t=1:size(tmat,2)
+ m=t-1;
+ if tmat(s,t,1)==1 % coeff requested
+ z2=LaguerreGaussianE([p,m,q,lambda],domain(:,:,1),domain(:,:,2),'pol');
+ coeffs(s,t,1)=overlap(z1,conj(z2),domain,domain(:,:,1),accuracy);
+ else
+ coeffs(s,t,1)=0;
+ end
+ if tmat(s,t,2)==1 % coeff requested
+ z2=LaguerreGaussianE([p,-m,q,lambda],domain(:,:,1),domain(:,:,2),'pol');
+ coeffs(s,t,2)=overlap(z1,conj(z2),domain,domain(:,:,1),accuracy);
+ else
+ coeffs(s,t,2)=0;
+ end
+ end
+ end
+end
\ No newline at end of file diff --git a/transverse/fitgaussianbeamrad.m b/transverse/fitgaussianbeamrad.m new file mode 100644 index 0000000..65e3cdb --- /dev/null +++ b/transverse/fitgaussianbeamrad.m @@ -0,0 +1,19 @@ +% Fits beamwidth versus position data to a Gaussian beam profile
+function f=fitgaussianbeamrad(z,w,lambda,varargin)
+
+if nargin>3, startparams=varargin{1};
+else
+ w0=1e-3;
+ z0=0;
+ startparams=[w0,z0];
+end
+if nargin>4, displayoptions=varargin{2}; else displayoptions=[1 0 0]; end
+
+f=genfit(z,w,startparams,@beamrad,'nonlin',displayoptions);
+
+
+function w=beamrad(params,z)
+w=beamradius([params,lambda],z); % Want 1/e^2 Power radius
+end
+
+end
\ No newline at end of file diff --git a/transverse/hermitepoly.m b/transverse/hermitepoly.m new file mode 100644 index 0000000..43b77d6 --- /dev/null +++ b/transverse/hermitepoly.m @@ -0,0 +1,18 @@ +%-----------------------------------
+% Hermite Polynomial
+% (See e.g. Arfken section 13.1)
+%
+% SYNTAX y=hermitepoly(n,x)
+%-----------------------------------
+
+function y=hermitepoly(n,x)
+
+m=[0:floor(n/2)];
+
+a=factorial(n-2*m);
+b=factorial(m);
+
+y=zeros(size(x));
+for s=1:length(m)
+ y = y + factorial(n) ./ a(s) ./ b(s) .* (-1).^m(s) * (2*x).^(n-2*m(s));
+end
diff --git a/transverse/overlap.m b/transverse/overlap.m new file mode 100644 index 0000000..f464303 --- /dev/null +++ b/transverse/overlap.m @@ -0,0 +1,83 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: overlap
+% AUTHOR: Andri M. Gretarsson
+% DATE: 7/10/04
+%
+% SYNTAX: coeff=overlap(z1,z2,domain <,metric,accuracy>)
+% <...> indicates optional arguments
+%
+% Calculates the overlap integral of the two functions specified.
+%
+% INPUT ARGUMENTS:
+% ----------------
+% z1 = The values of the first function. Can be a 1D vector or a 2D matrix
+% z2 = The values of the second function. z1 and z2 must be the same size.
+% domain = the domain values at which z1 and z2 are specified. If z1 and z2 are 1D
+% then domain is a nx2 matrix specifying one x,y pair for each value in z1
+% and z2. If z1 and z2 are 2D then domain is a nxmx2 array where
+% domain(:,:,1) and domain(:,:,2) are the x and y meshes corresponding
+% to the values in z1 and z2. These meshes are often generated using meshgrid.m.
+% metric = value, vector or matrix by which to multiply the elemental line or area dl or dS.
+% For example, in 2D polar coordinates the elemental area is dS=r*dr*dtheta so
+% metric should be specified as r.
+% accuracy = round results to the nearest increment of accuracy. For example, if
+% accuracy=0.3, then coeff 1.54 would be rounded to 1.5 while coeff=1.56 would be
+% rounded to 1.8.
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% coeff = the numerical result of the overlap integral.
+%
+% EXAMPLE 1 (cartesian):
+% [x,y]=meshgrid([0:0.01:2*pi],[0:0.01:2*pi]);
+% clear domain; domain(:,:,1)=x; domain(:,:,2)=y;
+% z1=sin(x+y); z2=sin(x+y);
+% coeff=overlap(z1,z2,domain,1,0.0001)
+% z1=sin(x+y); z2=cos(x+y);
+% coeff=overlap(z1,z2,domain,1,0.0001)
+%
+% EXAMPLE 2 (polar):
+% [r,theta]=meshgrid([0.01:0.01:1],[0:0.5:360]*pi/180);
+% clear domain; domain(:,:,1)=r; domain(:,:,2)=theta;
+% z1=r; z2=theta;
+% subplot(1,3,1); [x,y]=pol2cart(theta,r); h=pcolor(x,y,z1); set(h,'EdgeColor','none'); axis square;
+% subplot(1,3,2); [x,y]=pol2cart(theta,r); h=pcolor(x,y,z2); set(h,'EdgeColor','none'); axis square;
+% subplot(1,3,3); [x,y]=pol2cart(theta,r); h=pcolor(x,y,z1.*z2); set(h,'EdgeColor','none'); axis square;
+% coeff=overlap(z1,z2,domain,r,0.0001)
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+% SYNTAX: coeff=overlap(z1,z2,domain <,metric,accuracy>)
+%----------------------------------------------------------------------------------------------
+function coeff=overlap(z1,z2,domain,varargin)
+
+if nargin>=4;
+ metric=varargin{1};
+ if size(metric)==size(z1)
+ metric=metric(2:end,2:end);
+ end
+else
+ metric=1;
+end
+if nargin>=5; accuracy=varargin{2}; else accuracy=0; end
+if size(z1)~=size(z2), error('z1 and z2 must have same size'); end
+
+if min(size(z1))==1 & min(size(domain))==1 % 1D plot along axis
+ dl=metric.*(domain(2:end)-domain(1:end-1));
+ coeff=sum(z1(2:end).*z2(2:end).*dl)/2;
+end
+if min(size(z1))==1 & size(domain,2)==2 % 1D plot along arbitrary axis
+ dl=metric.*sqrt((domain(2:end,1)-domain(1:end-1,1)).^2 + (domain(2:end,2)-domain(1:end-1,2)).^2);
+ coeff=sum(z1(2:end).*z2(2:end).*dl)/2;
+end
+if min(size(z1))>= 2 % 2D plot over xy plane
+ coord1=domain(:,:,1); coord2=domain(:,:,2);
+ dcoord1=diff(coord1,1,2); dcoord1=dcoord1(2:end,:);
+ dcoord2=diff(coord2,1,1); dcoord2=dcoord2(:,2:end);
+ dS = metric.*dcoord1.*dcoord2;
+ coeff=sum(sum(z1(2:end,2:end).*z2(2:end,2:end).*dS));
+end
+
+if accuracy~=0
+ coeff=round(coeff/accuracy)*accuracy;
+end
\ No newline at end of file diff --git a/transverse/polarmesh.m b/transverse/polarmesh.m new file mode 100644 index 0000000..914c3e0 --- /dev/null +++ b/transverse/polarmesh.m @@ -0,0 +1,92 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: polarmesh
+% AUTHOR: Andri M. Gretarsson
+% DATE: 7/10/04
+%
+% SYNTAX: [rmesh,thetamesh,xmesh,ymesh]=polarmesh(r <,theta,rspacing>);
+% <...> indicates optional arguments
+%
+% Simplified way of generating arbitrary polar meshes.
+%
+% INPUT ARGUMENTS:
+% ----------------
+% r = specification for the r mesh. If r is a 1x3 vector the entries are assumed
+% to indicate r_start, r_end and n_rpts (number of r=const curves) and a set
+% of equally spaced or square spaced r=const mesh curves are generated. If r is
+% nx1, each entry is assumed to be a constant for which to generate an r=const
+% curve in the mesh. If r is 2D (nxn, where n>1) it is assumed to be already
+% in mesh form. In this case theta must be specified in mesh form and the
+% function does no more than pol2cart.
+% theta = specification for the theta mesh. If it is left blank the default theta is
+% used (a mesh with 1 line per theta from 0 to 360). If theta is 1x3, the
+% values are assumed to be [theta_start,theta_end,n_thetapts] and are used to
+% generate a set of n_thetapts equally spaced theta=const curves in the mesh.
+% If theta is [nx1] the one mesh curve is generated for each entry at the
+% specified angle.
+% rspacing = 'lin' for linear spacing of the r mesh
+% 'sqr' for r-line spacing proportional to r^2.
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% rmesh,thetamesh = r and theta meshes expressed in polar coordinates
+% xmesh,ymesh = the r and theta meshes expressed in cartesian coordinates.
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+% SYNTAX: [rmesh,thetamesh,xmesh,ymesh]=polarmesh(r <,theta,rspacing>);
+%----------------------------------------------------------------------------------------------
+function [rmesh,thetamesh,xmesh,ymesh]=polarmesh(r,varargin)
+
+if nargin>=2, theta=varargin{1}; defaulttheta=0; else defaulttheta=1; end
+if nargin>=3 && strcmp(varargin{2},'rsqr'), linflag=0; else linflag=1; end
+
+if min(size(r))==1 % r is 1D
+ if size(r,1)==1 && size(r,2)==3 % r=[r_start,r_end,rpts]
+ if linflag % linear r-spacing for mesh
+ rseed=(r(1):(r(2)-r(1))/(r(3)-1):r(2)).';
+ else % r-spacing proportional to r^2
+ rseed=(sqrt(r(1)):(sqrt(r(2))-sqrt(r(1)))/(r(3)-1):sqrt(r(2))).'.^2;
+ end
+ if defaulttheta==1 % default theta
+ thetaseed=(0:1:360).'*pi/180;
+ [thetamesh,rmesh]=meshgrid(thetaseed,rseed); % I like r to be arranged in the first position (r,theta)
+ thetamesh=transpose(thetamesh); % rather than (theta,r) so for consistency with (xmesh,ymesh)
+ rmesh=transpose(rmesh); % rmesh must have identical rows and thetamesh identical columns
+ else % not default theta
+ if size(theta,1)==1 && size(theta,2)==3 % theta =[theta_start,theta_end,theta_pts]
+ thetaseed=(theta(1):(theta(2)-theta(1))/(theta(3)-1):theta(2)).';
+ [thetamesh,rmesh]=meshgrid(thetaseed,rseed);
+ thetamesh=transpose(thetamesh);
+ rmesh=transpose(rmesh);
+ else % theta lines specified
+ [thetamesh,rmesh]=meshgrid(theta,rseed);
+ thetamesh=transpose(thetamesh);
+ rmesh=transpose(rmesh);
+ end
+ end
+ else % r lines specified
+ if size(r,1)<size(r,2), r=r.'; end % make r columnar
+ if defaulttheta==1 % default theta
+ thetaseed=(0:1:360).'*pi/180;
+ [thetamesh,rmesh]=meshgrid(thetaseed,r);
+ thetamesh=transpose(thetamesh);
+ rmesh=transpose(rmesh);
+ else % not default theta
+ if size(theta,1)==1 && size(theta,2)==3 % theta =[theta_start,theta_end,theta_pts]
+ thetaseed=(theta(1):(theta(2)-theta(1))/(theta(3)-1):theta(2)).';
+ [thetamesh,rmesh]=meshgrid(thetaseed,r);
+ thetamesh=transpose(thetamesh);
+ rmesh=transpose(rmesh);
+ else % theta lines specified
+ [thetamesh,rmesh]=meshgrid(theta,r); % theta and r must be the same length
+ thetamesh=transpose(thetamesh);
+ rmesh=transpose(rmesh);
+ end
+ end
+ end
+else
+ thetamesh = theta; % r and theta are already mesh grids
+ rmesh=r; % (and are the same size)
+end
+
+[xmesh,ymesh]=pol2cart(thetamesh,rmesh);
\ No newline at end of file diff --git a/transverse/recompose.m b/transverse/recompose.m new file mode 100644 index 0000000..0e23794 --- /dev/null +++ b/transverse/recompose.m @@ -0,0 +1,83 @@ +%----------------------------------------------------------------------------------------------
+% PROGRAM: recompose
+% AUTHOR: Andri M. Gretarsson
+% DATE: 7/10/04
+%
+% SYNTAX: z=recompose(domain,type,coeffs,[q <,lambda,accuracy>]);
+% <...> indicates optional arguments
+%
+% Calculates the complex field amplitude resulting from summing the specified terms of
+% a Hermite or Laguerre Gaussian mode expansion.
+%
+% INPUT ARGUMENTS:
+% ----------------
+% domain = the domain over which to do the recomposition. domain is a NxMx2 matrix
+% where domain(:,:,1) is the x mesh and domain(:,:,2) is the y mesh
+% (or r mesh and theta mesh respectively in the case of a Laguerre Gaussian
+% mode expansion).
+% type = 'hg' for a Hermite Gaussian mode expansion, and 'lg' for a Laguerre
+% Gaussian mode expansion.
+% coeffs = the matrix of coefficients in the same form as returned by decompose.m
+% q = the complex radius of curvature "q" of the Gaussian basis.
+% lambda = wavelength of the light in the Gaussian basis. Default is 1.064 microns
+% accuracy = only calculate the coefficients of the decomposition to th
+% specified accuracy. For example, if accuracy=0.3, then coeffs(i,j)= 1.4
+% would be rounded to 1.3.
+%
+% OUTPUT ARGUMENTS:
+% -----------------
+% z(i,j) = Resultant complex field of the recomposed modes at over the domain.
+%
+% Last updated: July 18, 2004 by AMG
+%----------------------------------------------------------------------------------------------
+% SYNTAX: z=recompose(domain,type,coeffs,[q <,lambda,accuracy>]);
+%----------------------------------------------------------------------------------------------
+
+function z=recompose(domain,type,coeffs,params)
+
+z=zeros(size(domain(:,:,1)));
+% HERMITE GAUSSIAN EXPANSION --------------------------------------------------------------------------
+if strcmpi(type,'hg')
+ q=params(1);
+ if length(params)>=2, lambda=params(2); else lambda=1.064e-6; end
+ if length(params)>=3, accuracy=params(3); else accuracy=1e-4; end
+ for s=1:size(coeffs,1)
+ l=s-1;
+ for t=1:size(coeffs,2)
+ m=t-1;
+ if abs(coeffs(s,t))>accuracy
+ z=z+HermiteGaussianE([l,m,q,lambda,coeffs(s,t)],domain(:,:,1),domain(:,:,2));
+ else
+ coeffs(s,t)=0;
+ end
+ end
+ end
+
+% LAGUERRE GAUSSIAN EXPANSION --------------------------------------------------------------------------
+elseif strcmpi(type,'lg')
+ q=params(1);
+ if length(params)>=2, lambda=params(2); else lambda=1.064e-6; end
+ if length(params)>=3, accuracy=params(3); else accuracy=1e-4; end
+ for s=1:size(coeffs,1)
+ p=s-1;
+ for t=1:size(coeffs,2)
+ m=t-1;
+ if abs(coeffs(s,t,1))>accuracy
+ z=z+LaguerreGaussianE([p,m,q,lambda,coeffs(s,t,1)],domain(:,:,1),domain(:,:,2));
+ else
+ coeffs(s,t,1)=0;
+ end
+ end
+ end
+ for s=1:size(coeffs,1)
+ p=s-1;
+ for t=2:size(coeffs,2) % skip terms with m=0
+ m=t-1;
+ if abs(coeffs(s,t,2))>accuracy
+ z=z+LaguerreGaussianE([p,-m,q,lambda,coeffs(s,t,2)],domain(:,:,1),domain(:,:,2));
+ else
+ coeffs(s,t,2)=0;
+ end
+ end
+ end
+end
\ No newline at end of file |