summaryrefslogtreecommitdiff
path: root/transverse/LaguerreGaussianE.m
blob: adc004a9dc97fd92505fa7dc4a46387bd9b32386 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
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))/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))/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