matlab中输入输入高斯光束的束腰与发散角及输出光束的光斑大小和发散角，求解衍射元件DOE的发散角，等参数的代码

以下是MATLAB代码示例：

% 输入参数
w0 = 1;         % 高斯光束束腰
theta0 = 0.1;   % 高斯光束发散角
wo = 10;        % 输出光束光斑大小
thetao = 0.2;   % 输出光束发散角

% 计算衍射元件DOE的发散角
lambda = 0.6328e-6;  % 波长
k = 2 * pi / lambda; % 波矢量
z = 1000;            % 衍射距离
L = 0.1;             % 衍射元件尺寸
M = 256;             % 衍射元件离散化点数
dx = L / M;          % 离散化步长
x = linspace(-L/2, L/2, M);
[X, Y] = meshgrid(x, x);
r = sqrt(X.^2 + Y.^2);
phi = atan2(Y, X);

Uin = w0^2 ./ (w0^2 + (r * tan(theta0)).^2) .* exp(-r.^2 ./ (w0^2 + (r * tan(theta0)).^2));  % 输入高斯光束
H = exp(1i * k * z * sqrt(1 - (lambda * r / (pi * wo^2)).^2)) ./ (1i * lambda * z * sqrt(1 - (lambda * r / (pi * wo^2)).^2)) .* exp(-1i * k * r.^2 / (2 * z)) .* exp(-1i * k * z) .* exp(1i * k * X.^2 / (2 * z)); % 衍射元件传递函数
Uout = H .* Uin;  % 衍射后的输出光束
Uout = Uout / max(max(abs(Uout))); % 归一化

% 计算输出光束发散角
Iout = abs(Uout).^2;  % 输出光强
Iout_norm = Iout / max(max(Iout));
x_int = trapz(Iout_norm, 1);  % 沿y轴积分
y_int = trapz(Iout_norm, 2);  % 沿x轴积分
x_int_norm = x_int / max(x_int);  % 归一化
y_int_norm = y_int / max(y_int);  % 归一化
x_int_norm_dB = 10 * log10(x_int_norm);
y_int_norm_dB = 10 * log10(y_int_norm);
x_int_max_dB = max(x_int_norm_dB);
y_int_max_dB = max(y_int_norm_dB);
theta_x = fwhm(x_int_norm_dB, x);  % 输出光束x方向发散角
theta_y = fwhm(y_int_norm_dB, x);  % 输出光束y方向发散角
theta = sqrt(theta_x^2 + theta_y^2);  % 输出光束总发散角

fprintf('衍射元件DOE的发散角为: %.3f\n', theta);

% 计算全宽度半最大值
function [w] = fwhm(y, x)
    half_max = max(y)/2;
    [~, idxs] = findpeaks(y);
    [~, idx1] = min(abs(y(1:idxs(1)) - half_max));
    [~, idx2] = min(abs(y(idxs(end):end) - half_max));
    idx2 = idx2 + idxs(end) - 1;
    w = abs(x(idx2) - x(idx1));
end

此代码用于计算高斯光束经过一段距离后经过某种衍射元件后的输出光束发散角。其中，输入高斯光束的束腰和发散角以及输出光束的光斑大小和发散角由用户输入，衍射元件尺寸和离散化点数固定为0.1和256。计算过程中使用了衍射元件传递函数，通过计算输出光束的强度分布沿x和y方向的积分值，进而求得输出光束的x和y方向的发散角，最终将两者平方和开方得到总发散角。