%----------------------------------------------------------------------------------------------
% 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(i*signm(u)*m(u).*theta).* LaguerrePoly([p(u),m(u)],2*r.^2/w(u).^2)...
.* exp( -i*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( -i*2*pi/lambda(u)*r.^2/2/q(u));
end
end
