涡旋光束阵列的产生【matlab代码】

  1. 求通过设计衍射光栅产生涡旋光束阵列的matlab程序,需要衍射级处相对光强的仿真图。

w0无法识别是怎么回事

clc
clear all
close all
%% 环形涡旋光束
N = 200;
lambda = 632e-9;    %波长为632nm
k = 2*pi/lambda;    %波数
w0 = 3;             %束腰半径
x = linspace(-10,10,N);
y = linspace(-10,10,N);
[X,Y] = meshgrid(x,y);
[theta,r] = cart2pol(X,Y);

beta = 50*pi/180;
figure;
for m = -4 : 4
    subplot(3,3,m+5)
    E1 = (r/w0).^abs(m).*exp(-r.^2/w0^2)*exp(1i*beta).*exp(-1i*m*theta);
    I1 = E1.*conj(E1);  I1 = I1/max(max(I1));
    h1 = pcolor(X,Y,I1);
    colorbar;
    set(h1,'edgecolor','none','facecolor','interp');
    title(['m = ',num2str(m)]);
    %colormap(gray);        %输出灰度图像
    axis square;
    %mesh(X,Y,I1)           %三维
end
suptitle('环形涡旋光束:不同拓扑荷数(m)')   %为图一添加总标题
%% 贝塞尔-高斯光束
N = 200;
lambda = 632e-9;    %波长为632nm
k = 2*pi/lambda;    %波数
w0 = 3;             %束腰半径
x = linspace(-2e-6,2e-6,N);
y = linspace(-2e-6,2e-6,N);
[X,Y] = meshgrid(x,y);
[theta,r] = cart2pol(X,Y);
figure;
alpha = sqrt(k^2-theta.^2);
for m = -4 : 4
    subplot(3,3,m+5)
    E2 = besselj(m,alpha.*r).*exp(-r.^2/w0^2).*exp(-1i*m*theta); %使用了matlab内置的贝塞尔函数
    I2 = E2.*conj(E2);  I2 = I2/max(max(I2));
    h2 = pcolor(X,Y,I2);
    colorbar;
    set(h2,'edgecolor','none','facecolor','interp');
    title(['m = ',num2str(m)]);
    %colormap(gray);        %输出灰度图像
    axis square;
    %mesh(X,Y,I2)           %三维
end
suptitle('贝塞尔-高斯光束:不同拓扑荷数(m)')   %为图二添加总标题
%% 拉盖尔-高斯光束
N = 200;
lambda = 632e-9;    %波长为632nm
k = 2*pi/lambda;    %波数
w0 = 3;          %光斑尺寸
x = linspace(-10,10,N);
y = linspace(-10,10,N);
[X,Y] = meshgrid(x,y);
[theta,r] = cart2pol(X,Y);

figure;
%p = 0, 1, 2;
p = 2;
for m = -4 : 4
    subplot(3,3,m+5)
   
    E3 = 1/w0*(sqrt(2)*r/w0).^abs(m).*laguerre(p,2*r.^2/w0^2).*exp(-r.^2/w0^2).*exp(1i*m*theta);
    I3 = E3.*conj(E3);  I3 = I3/max(max(I3));
    h3 = pcolor(X,Y,I3);
    colorbar;
    set(h3,'edgecolor','none','facecolor','interp');
    title(['m = ',num2str(m)]);
    %colormap(gray);        %输出灰度图像
    axis square;
    %mesh(X,Y,I3)           %三维
end
suptitle('拉盖尔-高斯光束:不同拓扑荷数(m)    p = 2')   %为图三添加总标题

%拉盖尔多项式
function result = laguerre(p,x)
result = 0;
switch p
    case 0
        result = 1;
    case 1
        for l = -p : p
            result = result+l+1-x;
        end
    case 2
        for l = -p : p
            result = result+0.5*(l+1)*(l+2)-(l+2)*x+0.5*x.^2;
        end
end
end