关于matlab的问题

Question

关于高斯光束整形的问题
需要把环形的高斯光束，加个偶级照明，截取两部分扇形，下面这个程序是高斯光束整形成环形的程序，我该在什么地方加个什么样的程序？最好是直接编写好的_(:з」∠)_最后成片是类似图中的样子，不过是有偏转，极角是45°

具体指标包括：入射光束为直径1.5mm平面波，光强为高斯分布，波长193nm；DOE通光口径为1.5mm，目标平面与DOE所在平面的距离为175mm，目标平面光场分布形式为偶极照明，外半径10.5mm，内半径7.5mm，极角45°，照明区域内要求光强均匀分布，光能利用率≥80%，光强不均匀性≤20%

%内能量超过90%的衍射圆
% 衍射距离。程序中的uint是毫米
% 定义变量
lamd=0.193*10^-3; % 入射波长
w=1.5; % 入射光束腰
R=1.5; % 元件的极限孔径
D=700; % 衍射距离
r1=1.5; % 圆的内半径
r2=3; % 圆的外半径
L0=10; % 水平长度（以毫米为单位）
k=2*pi/lamd; % 波数
N=256; % 采样点数
%%%%%%%%%%%%% 判断（分析变换）
% 判断是否满足菲涅尔变换条件
Judging=(sqrt(N*lamd*D)<=L0); % 计算一个非负数的平方根
if Judging==0
    disp('????ERROR');
    disp('......Fresnel Analytic Transformation is not satisfied'); % 菲涅尔变换不满足
    % break;
end
%%%%%%%%%%%%%
% 生成网格点
x11=linspace(-L0/2,L0/2,N); % 在[-L0/2, L0/2]范围内生成N个点
y11=linspace(-L0/2,L0/2,N); % 在[-L0/2, L0/2]范围内生成N个点
[x1,y1]=meshgrid(x11,y11); % 生成网格点
% 生成光强函数J1
J1=zeros(N);
for m=1:N
    for n=1:N
        if x1(m,n)^2+y1(m,n)^2<=R^2 % 如果点在圆内部，则取值为1，否则为0
            J1(m,n)=1;
        end
    end
end
% 生成干涉图的振幅函数A
A=exp(-(x1.^2+y1.^2)/w^2).*J1;
%%%%%
% 生成频率域上的网格点
fx=1/L0*(-N/2:N/2-1);
fy=1/L0*(-N/2:N/2-1);
[fx,fy]=meshgrid(fx,fy);
%%%%%%
% 生成干涉图的掩模函数JJ
JJ=zeros(N);
for m=1:N
    for n=1:N
        if x1(m,n)^2+y1(m,n)^2>=r1^2 && x1(m,n)^2+y1(m,n)^2<=r2^2 % 如果点在圆环内，则取值为1，否则为0
            JJ(m,n)=1;
        end
    end
end
a=sum(sum(A.^2))/sum(sum(JJ.^2)); % 计算放缩系数
J2=JJ*sqrt(a); % 对JJ进行缩放，使其平均光强等于A
%imagesc(J2);axis square;colormap(gray)图像SC（J2）；轴线正方形；颜色图（灰色）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% G-S
%pha0=load('C:\Documents and Settings\wjq\My Documents\MATLAB\pha90.txt');
pha0=2pirand(N)-pi; % 生成一个N×N大小的矩阵，矩阵中的每个元素都是在[-pi, pi]范围内随机生成的一个角度，作为初相。
M=0;CC=[];eta=0;
while eta<0.90 % 当eta小于0.9时，进行循环计算
M=M+1; % 循环次数自加1
U1=A.exp(ipha0); % 通过矩阵乘法生成初始的光强分布
temp1=(fft2(U1)); % 对初始光强分布进行二维傅里叶变换
temp2=exp(ikD*(1-lamd^2/2*(fx.^2+fy.^2))); % 计算平面波前传播的相位调制因子
U2=ifft2(temp1.temp2); % 对相位调制后的光强分布进行二维傅里叶逆变换，得到通过透镜后的光强分布
pha2=angle(U2); % 对通过透镜后的光强分布求相位，得到相位分布
A2=abs(U2); % 对通过透镜后的光强分布求幅值，得到幅值分布
I2=abs(U2).^2; % 对通过透镜后的光强分布求强度，得到强度分布
Err(M)=sum(sum((I2-J2.^2 ).^2)); % 计算输出光强与目标光强之间的误差
eta(M)=sum(sum(I2.JJ))/sum(sum(I2)); % 计算输出光强与目标光强之间的相似度
subplot(221);plotyy(1:M,Err,1:M,eta) % 绘制误差和相似度的变化曲线
title('Err (left) and eta (right)');xlabel('M');
subplot(222);imagesc(x11,y11,J2);axis square;title('Object function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') % 绘制目标光强分布的图像
subplot(223);imagesc(x11,y11,A2);axis square;title('Output function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') % 绘制透镜后的光强分布的图像
%%%%%%%%
U2=J2.exp(ipha2); % 将通过透镜后的光强分布的相位分布与目标光强分布相乘，得到下一次迭代的初始光强分布
tmep3=fft2(U2); % 对下一次迭代的初始光强分布进行二维傅里叶变换
temp4=exp(-ikD*(1-lamd^2/2*(fx.^2+fy.^2))); % 计算平面波后传播的相位调制因子
U11=ifft2(tmep3.*temp4); % 利用傅里叶变换将初始光强分布与相位调制因子相乘得到通过光阑后的光强分布，然后再进行傅里叶反变换得到下一次迭代的光强分布
subplot(224);imagesc(x11,y11,pha0);axis square; title('phase of DOE xlabel('x / mm');ylabel('y / mm') %绘图
pause(0.01)
end

Answer 1

我引用ChatGPT作答：
大概这样呢？不知道行不行。

lamd=0.193*10^-3; % 入射波长
w=1.5; % 入射光束腰
R=1.5; % 元件的极限孔径
D=175*10^-3; % 衍射距离
L0=10.5*10^-3; % 照明区域半径
r1=7.5*10^-3; % 内半径
r2=10.5*10^-3; % 外半径
%theta=45*pi/180; % 极角
theta = pi/4;
k=2*pi/lamd; % 波数
N=256;


% 生成坐标网格
x=linspace(-L0/2,L0/2,N);
y=linspace(-L0/2,L0/2,N);
[X,Y]=meshgrid(x,y);

% 生成圆形光阑
J1=zeros(N);
for m=1:N
    for n=1:N
        if X(m,n)^2+Y(m,n)^2<=R^2
            J1(m,n)=1;
        end
    end
end

% 生成高斯光束
A=exp(-(X.^2+Y.^2)/w^2).*J1;

% 生成偶极照明
J4=zeros(N);
for m=1:N
    for n=1:N
        if X(m,n)^2+Y(m,n)^2>=r1^2 && X(m,n)^2+Y(m,n)^2<=r2^2 && abs(atan2(Y(m,n),X(m,n)))<=theta/2
            J4(m,n)=1;
        end
    end
end

% 计算光阑和高斯光束的中心
J1_center = [round(N/2), round(N/2)];
A_center = [round(N/2), round(N/2)];

% 计算光阑和高斯光束需要平移的距离
J1_shift = J1_center - [find(J1(:,1),1,'first'), find(J1(1,:),1,'first')];
A_shift = A_center - [find(A(:,1),1,'first'), find(A(1,:),1,'first')];

% 平移光阑和高斯光束
J1 = circshift(J1, J1_shift);
A = circshift(A, A_shift);

% 生成扇形光阑
for m=1:N
    for n=1:N
        if J4(m,n) && abs(atan2(Y(m,n),X(m,n)))>theta/2
            J4(m,n)=0;
        end
    end
end

% 计算归一化系数
a=sum(sum(A.^2))/sum(sum(J4.^2));
J4=J4*sqrt(a);

% 变换到衍射面
A_fft=fftshift(fft2(A));
J4_fft=fftshift(fft2(J4));
x_fft=linspace(-1/(2*D),1/(2*D),N);
[X_fft,Y_fft]=meshgrid(x_fft,x_fft);
phase_term=exp(1i*k*(X_fft.^2+Y_fft.^2)*D);
B_fft=J4_fft.*phase_term.*A_fft;

% 变换回物面
B=ifft2(ifftshift(B_fft));

% 显示照明光斑
subplot(1,3,1);
imagesc(x,y,J4); axis square;
title('Illumination Pattern')
xlabel('x / mm'); ylabel('y / mm');


% 显示高斯光束
subplot(1,3,2);
imagesc(x11,y11,A);axis square;
title('Gaussian Beam')
xlabel('x / mm');ylabel('y / mm');
% 显示衍射图样
subplot(1,3,3);
imagesc(x,y,abs(B).^2); axis square;
title('Diffraction Pattern')
xlabel('x / mm'); ylabel('y / mm');

% 生成偶极照明和高斯光束的交叉项
J5 = zeros(N);
for m=1:N
for n=1:N
if ~isnan(J2(m,n)) && ~isnan(A(m,n)) && J2(m,n) && A(m,n)
J5(m,n) = 1;
end
end
end

% 显示交叉项
subplot(1,3,3);
imagesc(x11,y11,J5);axis square;
title('Cross Term')
xlabel('x / mm');ylabel('y / mm');