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

1 files changed, 197 insertions, 0 deletions

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