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
|
%----------------------------------------------------------------------------------------------
% PROGRAM: decompose
% AUTHOR: Andri M. Gretarsson, edited by Melissa Guidry
% DATE: 7/10/04
%
% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>],coordtype);
% <...> 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.
% coordtype = (string) label for the type of coordinates supplied in r and theta. If
% this argument 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.
%
% 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],'cart')
% 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>],coordtype);
%----------------------------------------------------------------------------------------------
function [coeffs,tmat]=decompose(z1,domain,type,terms,params,coordtype)
if strcmpi(coordtype,'cart')
metric = 1;
elseif strcmpi(coordtype, 'pol')
metric = domain(:,:,1);
end
% 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,metric,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),coordtype);
coeffs(s,t,1)=overlap(z1,conj(z2),domain,metric,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),coordtype);
coeffs(s,t,2)=overlap(z1,conj(z2),domain,metric,accuracy);
else
coeffs(s,t,2)=0;
end
end
end
end
|
