拉盖尔高斯光束_厄米高斯光束MATLAB仿真
激光原理by 贾而穑 130212114
厄米高斯光束MATLAB 仿真
其中主程序文件：plotHermiteGaussianBeams.m
子程序文件：HermitePoly.m
程序如下：
plotHermiteGaussianBeams.m
%-------------------------------------------------------------------------% % auther：Erse Jia
% Student ID 130212114
%-------------------------------------------------------------------------% %% Hermite Gaussian Beams
%% SET PARAMETERS
% Physical parameters
lambda = 500; % nm
k = 2*pi/lambda;
% The two parameters for the gaussian beam (and derived quantities) z0 = 1;
A0 = 1;
W0 = sqrt(lambda*z0/pi);
W = @(z) W0*sqrt(1+(z/z0)^2);
R = @(z) z*(1+(z/z0)^2);
Zeta = @(z) atan(z/z0);
% The coefficients for the Hermite-Gaussian (HG) beam of order (l,m)
A = [ 1 0 0 0;
1 1 0 0;
0 0 0 0;
0 0 .2 0];
% Display Parameters
res = 800;
z = 1e-9;
x = linspace(-2*W(z),2*W(z),res);
y = linspace(-2*W(z),2*W(z),res);
[X Y] = meshgrid(x,y);
X = X(:);
Y = Y(:);
%% RUN THE SIMULATION
% Preallocate Memory
U = zeros(length(X),1);
Utemp = zeros(length(X),1);
Utemp2 = zeros(length(X),1);
% Calculate Values that are independent of HG Polynomial order
lpf = exp(-1i*k*z - 1i*k*(X.^2 + Y.^2)/(2*R(z))); %lateral phase factor u = sqrt(2)*X/W(z);
v = sqrt(2)*Y/W(z);
for l = 1:size(A,1)
%if there are any terms of this order, calculate the x-HG (so you don't %need to repeat for each value of m
if sum(A(l,:) ~= 0) ~= 0
Utemp2 = (W0/W(z))*polyval(HermitePoly(l-1),u).*exp(-u.^2/2); else
continue;
end
for m = 1:size(A,2)
if A(l,m) ~= 0
Utemp = Utemp2.*(polyval(HermitePoly(m-1),v)).*exp(-v.^2/2); Utemp = A(l,m)*Utemp.*lpf*exp(1i*(l+m+1)*Zeta(z));
U = U + Utemp;
end
end
end
%% DRAW PLOTS
figure;
U = reshape(U,res,res);
imagesc(x,y,abs(U).^2);
axis square;
set(1,'color','w');
title('Hermite-Gaussian Beam of Order');
xlabel('x (nm)');
ylabel('y (nm)');
HermitePoly.m
%-------------------------------------------------------------------------% % HermitePoly.m by Erse Jia
% Student ID 130212114
% Given nonnegative integer n, compute the
% Hermite polynomial H_n. Return the result as a vector whose mth
% element is the coefficient of x^(n+1-m).
% polyval(HermitePoly(n),x) evaluates H_n(x).
%-------------------------------------------------------------------------%
function hk = HermitePoly(n)
if n==0
hk = 1;
elseif n==1
hk = [2 0];
else
hkm2 = zeros(1,n+1);
hkm2(n+1) = 1;
hkm1 = zeros(1,n+1);
hkm1(n) = 2;
for k=2:n
hk = zeros(1,n+1);
for e=n-k+1:2:n
hk(e) = 2*(hkm1(e+1) - (k-1)*hkm2(e));
end
hk(n+1) = -2*(k-1)*hkm2(n+1);
if k
hkm2 = hkm1;
hkm1 = hk;
end
end
end
结果：
拉盖尔高斯光束MATLAB 仿真
主程序文件：DrawtheLaguerreGaussbeam.m
子程序文件：LG.m
DrawtheLaguerreGaussbeam.m
%-------------------------------------------------------------------------% % auther：Erse Jia
% Student ID 130212114
%-------------------------------------------------------------------------% clear all;
close all; clc
params = [0 0 1];
% Use function handle
u0 = @(rho, phi)LG(params, rho, phi);
R = @(x, y)(x=0);
u = @(rho, phi)R(rho.*cos(phi), rho.*sin(phi)).*u0(rho, phi);
[X, Y] = meshgrid(linspace(-5, 5, 200));
Rho = sqrt(X.^2 + Y.^2);
Phi = atan(Y./X);
figure(1)
set(1,'color','w');
Z = u(Rho, Phi);
surf(X, Y, -Z)
shading interp
set(gca,'box','on');
grid off;
xlabel('x position');
ylabel('y position');
zlabel('z');
colorbar;
LG.m
function y = LG( params, rho, phi )
m = abs(params(1));
p = params(2);
w = params(3);
if w==0
msgbox('params(0) can not be equal to 0');
end
t = rho./w;
y = sqrt(2*factorial(p)/pi/factorial(m+p))/w.* (sqrt(2).*t).^m ...
.* L([p m], 2*t.^2).* exp(-t.^2 + 1i*m*phi);
function y = L(params, x)
fact = @(x)arrayfun(@factorial, x);
n = params(1); % p
k = params(2); % m
m = 0:n;
a = factorial(n+k)*ones(1,length(m));
b = fact(n-m);
c = fact(k+m);
d = fact(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
end
end
结果：