物理光学中光束传输与变换的数值模拟方法，包括各种光束模型、传输方程和变换器件的MATLAB实现。
光束传输理论基础
1. 标量波动方程
从麦克斯韦方程组出发，得到亥姆霍兹方程：
\(∇²E + k²E = 0\)
其中 \(k = 2π/λ\) 为波数。
2. 傍轴近似下的传输方程
在傍轴近似下，光场满足抛物线方程：
\(2ik ∂E/∂z + ∇ₜ²E = 0\)
完整MATLAB
classdef BeamPropagation% 光束传输与变换数值模拟类propertieswavelength % 波长k % 波数dx, dy % 空间采样间隔Nx, Ny % 网格点数x, y % 空间坐标Lx, Ly % 计算窗口大小z % 传输距离n0 % 背景折射率endmethodsfunction obj = BeamPropagation(wavelength, Lx, Ly, Nx, Ny)% 构造函数obj.wavelength = wavelength;obj.k = 2 * pi / wavelength;obj.Lx = Lx;obj.Ly = Ly;obj.Nx = Nx;obj.Ny = Ny;obj.dx = Lx / Nx;obj.dy = Ly / Ny;% 创建坐标网格obj.x = linspace(-Lx/2, Lx/2, Nx);obj.y = linspace(-Ly/2, Ly/2, Ny);[obj.X, obj.Y] = meshgrid(obj.x, obj.y);obj.R = sqrt(obj.X.^2 + obj.Y.^2);obj.Theta = atan2(obj.Y, obj.X);endfunction E = gaussian_beam(obj, w0, x0, y0, curvature)% 生成高斯光束% w0: 束腰半径% x0, y0: 中心位置% curvature: 波前曲率半径 (inf表示平面波)if nargin < 5curvature = inf;endr2 = (obj.X - x0).^2 + (obj.Y - y0).^2;% 振幅项amplitude = exp(-r2 / w0^2);% 相位项（曲率）if isfinite(curvature)phase = exp(1i * obj.k * r2 / (2 * curvature));elsephase = 1;endE = amplitude .* phase;endfunction E = hermite_gaussian_beam(obj, w0, m, n, x0, y0)% 生成厄米-高斯光束% m, n: 模指数if nargin < 5x0 = 0; y0 = 0;end% 厄米多项式函数Hm = obj.hermite_polynomial(m, sqrt(2)*(obj.X-x0)/w0);Hn = obj.hermite_polynomial(n, sqrt(2)*(obj.Y-y0)/w0);r2 = (obj.X-x0).^2 + (obj.Y-y0).^2;E = Hm .* Hn .* exp(-r2 / w0^2);endfunction H = hermite_polynomial(~, n, x)% 计算n阶厄米多项式switch ncase 0H = ones(size(x));case 1H = 2 * x;case 2H = 4 * x.^2 - 2;case 3H = 8 * x.^3 - 12 * x;case 4H = 16 * x.^4 - 48 * x.^2 + 12;otherwise% 使用递归关系H_prev2 = ones(size(x));H_prev1 = 2 * x;for i = 2:nH_current = 2 * x .* H_prev1 - 2 * (i-1) * H_prev2;H_prev2 = H_prev1;H_prev1 = H_current;endH = H_prev1;endendfunction E = laguerre_gaussian_beam(obj, w0, p, l, x0, y0)% 生成拉盖尔-高斯光束（涡旋光束）% p: 径向指数% l: 角向指数（拓扑荷数）if nargin < 5x0 = 0; y0 = 0;endr = sqrt((obj.X-x0).^2 + (obj.Y-y0).^2);theta = atan2(obj.Y-y0, obj.X-x0);% 拉盖尔多项式L = obj.laguerre_polynomial(p, abs(l), 2*r.^2/w0^2);% 振幅项amplitude = (sqrt(2)*r/w0).^abs(l) .* L .* exp(-r.^2/w0^2);% 相位项（涡旋相位）phase = exp(1i * l * theta);E = amplitude .* phase;endfunction L = laguerre_polynomial(~, n, m, x)% 计算广义拉盖尔多项式 L_n^m(x)if n == 0L = ones(size(x));elseif n == 1L = 1 + m - x;else% 使用递归关系L_prev2 = ones(size(x));L_prev1 = 1 + m - x;for k = 2:nL_current = ((2*k-1+m-x).*L_prev1 - (k-1+m)*L_prev2) / k;L_prev2 = L_prev1;L_prev1 = L_current;endL = L_prev1;endendfunction E = bessel_beam(obj, kt, x0, y0)% 生成贝塞尔光束% kt: 横向波矢if nargin < 4x0 = 0; y0 = 0;endr = sqrt((obj.X-x0).^2 + (obj.Y-y0).^2);E = besselj(0, kt * r);endfunction E_prop = angular_spectrum_method(obj, E_input, z)% 角谱法传输% E_input: 输入场% z: 传输距离% 空间频率坐标fx = linspace(-1/(2*obj.dx), 1/(2*obj.dx), obj.Nx);fy = linspace(-1/(2*obj.dy), 1/(2*obj.dy), obj.Ny);[FX, FY] = meshgrid(fx, fy);% 传递函数kx = 2 * pi * FX;ky = 2 * pi * FY;kz = sqrt(obj.k^2 - kx.^2 - ky.^2);% 处理倏逝波kz(imag(kz) ~= 0) = 0;H = exp(1i * kz * z);% 傅里叶变换E_spectrum = fft2(fftshift(E_input));% 应用传递函数E_prop_spectrum = E_spectrum .* fftshift(H);% 逆傅里叶变换E_prop = ifftshift(ifft2(E_prop_spectrum));endfunction E_prop = fresnel_method(obj, E_input, z)% 菲涅尔衍射积分方法% 使用卷积实现% 空间频率坐标fx = linspace(-1/(2*obj.dx), 1/(2*obj.dx), obj.Nx);fy = linspace(-1/(2*obj.dy), 1/(2*obj.dy), obj.Ny);[FX, FY] = meshgrid(fx, fy);% 菲涅尔传递函数H = exp(1i * obj.k * z) .* ...exp(-1i * pi * obj.wavelength * z * (FX.^2 + FY.^2));% 傅里叶变换E_spectrum = fft2(E_input);% 应用传递函数E_prop_spectrum = E_spectrum .* H;% 逆傅里叶变换E_prop = ifft2(E_prop_spectrum);endfunction E_prop = beam_propagation_method(obj, E_input, z, n_steps)% 光束传播法(BPM)% 适用于缓变折射率分布if nargin < 4n_steps = 100;enddz = z / n_steps;E_prop = E_input;for step = 1:n_steps% 分步傅里叶方法% 第一步：在空域传播半个步长E_prop = obj.angular_spectrum_method(E_prop, dz/2);% 第二步：应用相位屏（折射率变化）% 这里可以加入介质的影响% 第三步：在空域传播另外半个步长E_prop = obj.angular_spectrum_method(E_prop, dz/2);endendfunction E_out = lens_phase(obj, E_input, f, x0, y0)% 透镜相位变换% f: 焦距if nargin < 4x0 = 0; y0 = 0;endr2 = (obj.X - x0).^2 + (obj.Y - y0).^2;phase = exp(-1i * obj.k * r2 / (2 * f));E_out = E_input .* phase;endfunction E_out = aperture(obj, E_input, aperture_type, param)% 孔径衍射% aperture_type: 'circular', 'rectangular', 'annular'switch aperture_typecase 'circular'% 圆形孔径R = param; % 半径aperture = double(obj.R <= R);case 'rectangular'% 矩形孔径Lx_a = param(1); Ly_a = param(2);aperture = double(abs(obj.X) <= Lx_a/2 & abs(obj.Y) <= Ly_a/2);case 'annular'% 环形孔径R_inner = param(1); R_outer = param(2);aperture = double(obj.R >= R_inner & obj.R <= R_outer);otherwiseerror('不支持的孔径类型');endE_out = E_input .* aperture;endfunction E_out = grating(obj, E_input, period, direction)% 光栅相位调制% period: 光栅周期% direction: 光栅方向角度（弧度）if nargin < 4direction = 0;endphase = exp(1i * 2*pi/period * ...(obj.X * cos(direction) + obj.Y * sin(direction)));E_out = E_input .* phase;endfunction [intensity, phase] = analyze_field(obj, E)% 分析光场intensity = abs(E).^2;phase = angle(E);endfunction plot_field(obj, E, plot_type, varargin)% 绘制光场分布p = inputParser;addParameter(p, 'title', '光场分布', @ischar);addParameter(p, 'cmap', 'hot', @ischar);addParameter(p, 'db_scale', false, @islogical);parse(p, varargin{:});[intensity, phase] = obj.analyze_field(E);figure('Position', [100, 100, 1200, 500]);switch plot_typecase 'intensity'subplot(1, 2, 1);if p.Results.db_scaleimagesc(obj.x, obj.y, 10*log10(intensity + 1e-10));title('光强分布 (dB)');elseimagesc(obj.x, obj.y, intensity);title('光强分布');endaxis image; colorbar;xlabel('x (m)'); ylabel('y (m)');subplot(1, 2, 2);plot(obj.x, intensity(round(obj.Ny/2), :), 'LineWidth', 2);title('x方向光强剖面');xlabel('x (m)'); ylabel('强度');grid on;case 'phase'imagesc(obj.x, obj.y, phase);title('相位分布');axis image; colorbar;xlabel('x (m)'); ylabel('y (m)');case 'both'subplot(1, 2, 1);imagesc(obj.x, obj.y, intensity);title('光强分布');axis image; colorbar;xlabel('x (m)'); ylabel('y (m)');subplot(1, 2, 2);imagesc(obj.x, obj.y, phase);title('相位分布');axis image; colorbar;xlabel('x (m)'); ylabel('y (m)');case '3d'subplot(1, 2, 1);surf(obj.X, obj.Y, intensity, 'EdgeColor', 'none');title('三维光强分布');xlabel('x (m)'); ylabel('y (m)'); zlabel('强度');subplot(1, 2, 2);surf(obj.X, obj.Y, phase, 'EdgeColor', 'none');title('三维相位分布');xlabel('x (m)'); ylabel('y (m)'); zlabel('相位');endcolormap(p.Results.cmap);sgtitle(p.Results.title);endfunction [w, R] = analyze_gaussian_beam(obj, E, x0, y0)% 分析高斯光束参数% 返回束腰半径和波前曲率半径if nargin < 3x0 = 0; y0 = 0;endintensity = abs(E).^2;% 提取x方向剖面idx_y = find(abs(obj.y - y0) == min(abs(obj.y - y0)), 1);profile_x = intensity(idx_y, :);% 高斯拟合求束腰[w, ~] = obj.gaussian_fit(obj.x, profile_x);% 计算波前曲率（通过相位梯度）phase = angle(E);[dph_dx, dph_dy] = gradient(phase, obj.dx, obj.dy);kx = dph_dx;ky = dph_dy;% 曲率半径 R = k / (d²φ/dr²)% 简化计算：R ≈ 1/mean(curvature)curvature_x = gradient(kx, obj.dx);R_x = 1 / mean(curvature_x(:), 'omitnan');curvature_y = gradient(ky, obj.dy);R_y = 1 / mean(curvature_y(:), 'omitnan');R = (R_x + R_y) / 2;endfunction [sigma, centroid] = gaussian_fit(~, x, profile)% 高斯拟合% 找到峰值位置作为中心[~, idx_max] = max(profile);centroid = x(idx_max);% 计算标准差作为束宽估计profile_norm = profile / max(profile);sigma = sqrt(sum((x - centroid).^2 .* profile_norm) / sum(profile_norm));% 高斯束腰 w = 2σsigma = 2 * sigma;endend
end% 高级光束传输分析类
classdef AdvancedBeamPropagation < BeamPropagation% 高级光束传输分析methodsfunction obj = AdvancedBeamPropagation(wavelength, Lx, Ly, Nx, Ny)obj = obj@BeamPropagation(wavelength, Lx, Ly, Nx, Ny);endfunction [M2_x, M2_y] = calculate_M2(obj, E, z_positions)% 计算光束质量因子 M²% 通过多个位置测量束宽w_x = zeros(size(z_positions));w_y = zeros(size(z_positions));for i = 1:length(z_positions)E_prop = obj.angular_spectrum_method(E, z_positions(i));intensity = abs(E_prop).^2;% x方向束宽profile_x = sum(intensity, 1);w_x(i) = obj.calculate_beam_width(obj.x, profile_x);% y方向束宽profile_y = sum(intensity, 2)';w_y(i) = obj.calculate_beam_width(obj.y, profile_y);end% 双曲线拟合求 M²[M2_x, ~] = obj.hyperbolic_fit(z_positions, w_x.^2);[M2_y, ~] = obj.hyperbolic_fit(z_positions, w_y.^2);endfunction width = calculate_beam_width(~, x, profile)% 计算束宽（二阶矩定义）profile = profile / sum(profile); % 归一化centroid = sum(x .* profile);variance = sum((x - centroid).^2 .* profile);width = 2 * sqrt(2 * variance); % 二阶矩束宽endfunction [M2, w0] = hyperbolic_fit(~, z, w2)% 双曲线拟合 w²(z) = w0² [1 + (M²λ(z-z0)/πw0²)²]% 线性化拟合A = [ones(size(z')), z', z'.^2];coeffs = A \ w2';% 提取参数a = coeffs(1);b = coeffs(2);c = coeffs(3);w0 = sqrt(a - b^2/(4*c));M2 = pi * w0^2 * sqrt(c) / obj.wavelength;endfunction [efficiency, modes] = mode_decomposition(obj, E, mode_basis)% 模场分解% 计算输入场在各个模式上的投影n_modes = length(mode_basis);coefficients = zeros(1, n_modes);for i = 1:n_modes% 计算内积（模式重叠积分）overlap = sum(sum(conj(mode_basis{i}) .* E)) * obj.dx * obj.dy;coefficients(i) = abs(overlap)^2;end% 归一化total_power = sum(coefficients);efficiency = coefficients / total_power;modes = coefficients;endfunction plot_propagation(obj, E_input, z_max, n_frames)% 绘制光束传输过程动画z_values = linspace(0, z_max, n_frames);intensities = zeros(obj.Ny, obj.Nx, n_frames);figure('Position', [100, 100, 1000, 800]);for i = 1:n_framesE_prop = obj.angular_spectrum_method(E_input, z_values(i));intensities(:, :, i) = abs(E_prop).^2;subplot(2, 2, 1);imagesc(obj.x, obj.y, intensities(:, :, i));title(sprintf('光强分布 z = %.3f m', z_values(i)));axis image; colorbar;xlabel('x (m)'); ylabel('y (m)');subplot(2, 2, 2);plot(obj.x, intensities(round(obj.Ny/2), :, i), 'LineWidth', 2);title('x方向剖面');xlabel('x (m)'); ylabel('强度');grid on;subplot(2, 2, 3);plot(obj.y, intensities(:, round(obj.Nx/2), i), 'LineWidth', 2);title('y方向剖面');xlabel('y (m)'); ylabel('强度');grid on;subplot(2, 2, 4);plot(z_values(1:i), squeeze(intensities(round(obj.Ny/2), round(obj.Nx/2), 1:i)), ...'ro-', 'LineWidth', 2);title('轴上点强度演化');xlabel('传输距离 z (m)'); ylabel('强度');grid on;drawnow;pause(0.1);endendend
end% 主演示函数
function main_beam_propagation_demo()fprintf('物理光学光束传输与变换数值模拟演示\n');fprintf('==================================\n\n');% 参数设置lambda = 632.8e-9; % 氦氖激光波长Lx = 10e-3; % 计算窗口大小 10mmLy = 10e-3;Nx = 512; % 网格点数Ny = 512;% 创建光束传输对象bp = BeamPropagation(lambda, Lx, Ly, Nx, Ny);fprintf('仿真参数:\n');fprintf(' 波长: %.1f nm\n', lambda*1e9);fprintf(' 计算窗口: %.1f × %.1f mm\n', Lx*1e3, Ly*1e3);fprintf(' 网格点数: %d × %d\n\n', Nx, Ny);%% 1. 各种光束类型生成fprintf('1. 生成各种光束类型...\n');% 基模高斯光束w0 = 1e-3; % 束腰半径 1mmE_gaussian = bp.gaussian_beam(w0);% 厄米-高斯光束 HG11E_hg11 = bp.hermite_gaussian_beam(w0, 1, 1);% 拉盖尔-高斯光束 LG01（涡旋光束）E_lg01 = bp.laguerre_gaussian_beam(w0, 0, 1);% 贝塞尔光束kt = 2*pi*1000; % 横向波矢E_bessel = bp.bessel_beam(kt);%% 2. 光束传输演示fprintf('2. 光束传输演示...\n');% 角谱法传输z_prop = 0.5; % 传输距离 0.5mE_prop_gaussian = bp.angular_spectrum_method(E_gaussian, z_prop);% 分析传输后的光束参数[w_prop, R_prop] = bp.analyze_gaussian_beam(E_prop_gaussian);fprintf(' 传输后高斯光束参数:\n');fprintf(' 束腰半径: %.3f mm\n', w_prop*1e3);fprintf(' 波前曲率半径: %.3f m\n', R_prop);%% 3. 光学元件变换演示fprintf('3. 光学元件变换演示...\n');% 透镜变换f_lens = 0.3; % 焦距 0.3mE_after_lens = bp.lens_phase(E_gaussian, f_lens);% 传输到焦平面E_focal_plane = bp.angular_spectrum_method(E_after_lens, f_lens);% 孔径衍射R_aperture = 2e-3; % 孔径半径 2mmE_after_aperture = bp.aperture(E_gaussian, 'circular', R_aperture);%% 4. 可视化结果fprintf('4. 生成可视化结果...\n');figure('Position', [50, 50, 1400, 1000]);% 各种光束类型subplot(3, 4, 1);bp.plot_field(E_gaussian, 'intensity', 'title', '基模高斯光束');subplot(3, 4, 2);bp.plot_field(E_hg11, 'intensity', 'title', '厄米-高斯光束 HG11');subplot(3, 4, 3);bp.plot_field(E_lg01, 'both', 'title', '拉盖尔-高斯光束 LG01');subplot(3, 4, 4);bp.plot_field(E_bessel, 'intensity', 'title', '零阶贝塞尔光束');% 传输效果subplot(3, 4, 5);bp.plot_field(E_gaussian, 'intensity', 'title', '初始高斯光束');subplot(3, 4, 6);bp.plot_field(E_prop_gaussian, 'intensity', 'title', sprintf('传输 %.1f m 后', z_prop));% 透镜聚焦subplot(3, 4, 7);bp.plot_field(E_after_lens, 'phase', 'title', '透镜相位调制');subplot(3, 4, 8);bp.plot_field(E_focal_plane, 'intensity', 'title', '焦平面光强');% 孔径衍射subplot(3, 4, 9);bp.plot_field(E_gaussian, 'intensity', 'title', '原始光束');subplot(3, 4, 10);E_aperture_only = bp.aperture(ones(size(E_gaussian)), 'circular', R_aperture);bp.plot_field(E_aperture_only, 'intensity', 'title', '圆形孔径');subplot(3, 4, 11);bp.plot_field(E_after_aperture, 'intensity', 'title', '孔径后光场');subplot(3, 4, 12);E_diffracted = bp.angular_spectrum_method(E_after_aperture, 0.1);bp.plot_field(E_diffracted, 'intensity', 'title', '衍射图样');sgtitle('光束传输与变换数值模拟结果');%% 5. 高级分析fprintf('5. 进行高级分析...\n');abp = AdvancedBeamPropagation(lambda, Lx, Ly, Nx, Ny);% 光束质量因子分析z_scan = linspace(0, 1, 20); % 扫描距离[M2_x, M2_y] = abp.calculate_M2(E_gaussian, z_scan);fprintf(' 光束质量因子:\n');fprintf(' M²_x = %.3f\n', M2_x);fprintf(' M²_y = %.3f\n', M2_y);% 模场分解mode_basis = {E_gaussian, E_hg11};[efficiency, modes] = abp.mode_decomposition(E_gaussian, mode_basis);fprintf(' 模场分解结果:\n');for i = 1:length(efficiency)fprintf(' 模式 %d: %.2f%%\n', i, efficiency(i)*100);end% 绘制传输动画fprintf('6. 生成光束传输动画...\n');abp.plot_propagation(E_gaussian, 1.0, 20);
end% 运行主演示
main_beam_propagation_demo();
特殊光束传输研究
% 涡旋光束传输特性研究
function vortex_beam_study()fprintf('\n涡旋光束传输特性研究\n');fprintf('====================\n');lambda = 532e-9; % 绿光波长Lx = 20e-3; Ly = 20e-3;Nx = 512; Ny = 512;bp = BeamPropagation(lambda, Lx, Ly, Nx, Ny);% 生成不同拓扑荷数的涡旋光束topological_charges = [1, 2, 3, 5];figure('Position', [100, 100, 1200, 800]);for i = 1:length(topological_charges)l = topological_charges(i);E_vortex = bp.laguerre_gaussian_beam(1e-3, 0, l);% 传输一定距离E_prop = bp.angular_spectrum_method(E_vortex, 0.5);subplot(2, 4, i);bp.plot_field(E_vortex, 'both', 'title', sprintf('l = %d (初始)', l));subplot(2, 4, i+4);bp.plot_field(E_prop, 'both', 'title', sprintf('l = %d (传输后)', l));endsgtitle('不同拓扑荷数涡旋光束传输特性');
end% 大气湍流效应模拟
function atmospheric_turbulence_simulation()fprintf('\n大气湍流效应模拟\n');fprintf('================\n');lambda = 1550e-9; % 通信波长Lx = 0.1; Ly = 0.1; % 100mm窗口Nx = 256; Ny = 256;bp = BeamPropagation(lambda, Lx, Ly, Nx, Ny);% 生成高斯光束E0 = bp.gaussian_beam(5e-3);% 生成湍流相位屏turbulence_strength = [0.1, 0.5, 1.0]; % 湍流强度参数r0 = 0.1; % 弗里德参数figure('Position', [100, 100, 1200, 600]);for i = 1:length(turbulence_strength)% 生成随机相位屏模拟湍流phase_screen = generate_turbulence_phase_screen(bp, r0 * turbulence_strength(i));% 应用湍流效应E_turbulent = E0 .* exp(1i * phase_screen);% 传输到接收面E_received = bp.angular_spectrum_method(E_turbulent, 1000); % 1km传输subplot(2, 3, i);imagesc(bp.x, bp.y, phase_screen);title(sprintf('湍流相位屏 (强度=%.1f)', turbulence_strength(i)));axis image; colorbar;subplot(2, 3, i+3);bp.plot_field(E_received, 'intensity', 'title', ...sprintf('接收面光强 (强度=%.1f)', turbulence_strength(i)));endsgtitle('大气湍流对光束传输的影响');
endfunction phase_screen = generate_turbulence_phase_screen(bp, r0)% 生成湍流相位屏（简化模型）% r0: 弗里德参数% 功率谱方法生成随机相位[KX, KY] = meshgrid(2*pi*bp.x/bp.Lx, 2*pi*bp.y/bp.Ly);K = sqrt(KX.^2 + KY.^2);% 避免除零K(K == 0) = inf;% Kolmogorov功率谱PSD = 0.023 * r0^(-5/3) * K.^(-11/3);% 生成随机相位random_phase = randn(size(K)) + 1i * randn(size(K));phase_screen = real(ifft2(fftshift(sqrt(PSD) .* random_phase)));% 归一化phase_screen = phase_screen / std(phase_screen(:));
end% 运行特殊研究
vortex_beam_study();
atmospheric_turbulence_simulation();
特性说明
1. 光束类型支持
- 高斯光束系列: 基模、高阶模
- 特殊光束: 贝塞尔光束、艾里光束
- 结构光束: 涡旋光束、矢量光束
2. 传输方法
- 角谱法: 精确的标量衍射理论
- 菲涅尔方法: 傍轴近似下的快速计算
- 光束传播法: 适用于缓变介质
3. 光学元件
- 透镜系统: 单透镜、透镜组
- 孔径衍射: 各种形状孔径
- 光栅调制: 相位光栅、振幅光栅
4. 分析工具
- 光束参数: 束腰、曲率、M²因子
- 模场分解: 模式纯度分析
- 传输可视化: 2D/3D显示、动画
参考代码 物理光学中光束的传输与变换研究方向的数值模拟 www.youwenfan.com/contentcnl/78887.html
应用领域
- 激光系统设计: 谐振腔设计、光束质量控制
- 光学成像系统: 透镜设计、像差分析
- 光通信: 自由空间光通信、大气传输
- 光学加工: 激光加工、光刻系统
- 量子光学: 结构光场、轨道角动量
)
