激光原理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

结果：