MATLAB:从振幅到能流——菲涅尔公式的功率传输可视化
1. 菲涅尔公式与能量传输的物理本质当一束光从空气照射到玻璃表面时你会发现部分光线被反射部分光线进入玻璃继续传播。这个看似简单的现象背后隐藏着法国物理学家奥古斯丁·菲涅尔在19世纪发现的深刻规律。菲涅尔公式不仅能够精确计算反射和透射的振幅比例更重要的是揭示了能量在光学界面处的分配机制。在实际应用中我们往往更关心的是能量而非振幅。因为光的强度单位面积上的功率与电场振幅的平方成正比。这就引出了能流反射率R和透射率T的概念R|r|²表示反射能量占总入射能量的比例T(n₂cosθₜ/n₁cosθᵢ)·|t|²则表示透射能量的比例。这两个量才是真正决定有多少光被反射回来多少光继续前进的关键指标。我记得第一次用MATLAB绘制能流曲线时发现当入射角达到某个特定值布儒斯特角时p偏振光的反射率竟然降为零。这意味着所有能量都进入了第二种介质这个反直觉的现象让我对菲涅尔公式的精确性产生了深刻印象。2. MATLAB实现光疏到光密介质的能流计算让我们以空气n₁1.0到玻璃n₂1.5的界面为例看看如何用MATLAB实现能流计算。首先需要明确几个关键公式s偏振光的反射系数rₛ (n₁cosθᵢ - n₂cosθₜ)/(n₁cosθᵢ n₂cosθₜ)p偏振光的反射系数rₚ (n₂cosθᵢ - n₁cosθₜ)/(n₂cosθᵢ n₁cosθₜ)透射系数则根据能量守恒由t√(1-r²)得出% 光疏到光密介质能流计算 n1 1.0; n2 1.5; % 折射率 theta_i 0:0.1:90; % 入射角范围 a theta_i*pi/180; % 角度转弧度 % 计算折射角需要处理全反射情况 theta_t asin(n1/n2*sin(a)); % s偏振反射系数 rs (n1*cos(a) - n2*cos(theta_t))./(n1*cos(a) n2*cos(theta_t)); % p偏振反射系数 rp (n2*cos(a) - n1*cos(theta_t))./(n2*cos(a) n1*cos(theta_t)); % 能流反射率 Rs abs(rs).^2; Rp abs(rp).^2; % 能流透射率考虑折射率修正因子 Ts (n2*cos(theta_t)./(n1*cos(a))).*(1 - Rs); Tp (n2*cos(theta_t)./(n1*cos(a))).*(1 - Rp);绘制出的曲线会显示几个关键特征在小角度入射时反射率较低随着角度增大反射率逐渐升高对于p偏振光在布儒斯特角约56.3°处反射率降为零。特别值得注意的是能流的RT1验证了能量守恒。3. 光密到光疏介质的全反射现象当光从玻璃n₁1.5射向空气n₂1.0时情况变得更有趣。此时存在一个临界角θ_carcsin(n₂/n₁)≈41.8°当入射角超过这个值时会发生全反射现象。% 光密到光疏介质 n1 1.5; n2 1.0; theta_i 0:0.1:90; a theta_i*pi/180; % 处理全反射情况 sin_theta_t n1/n2*sin(a); theta_t real(asin(sin_theta_t)); % 取实部避免复数 % 反射系数计算 rs (n1*cos(a) - n2*cos(theta_t))./(n1*cos(a) n2*cos(theta_t)); rp (n2*cos(a) - n1*cos(theta_t))./(n2*cos(a) n1*cos(theta_t)); % 能流计算 Rs abs(rs).^2; Rp abs(rp).^2; Ts (n2*cos(theta_t)./(n1*cos(a))).*(1 - Rs); Tp (n2*cos(theta_t)./(n1*cos(a))).*(1 - Rp); % 修正全反射区域的透射率 critical_angle asind(n2/n1); Ts(theta_i critical_angle) 0; Tp(theta_i critical_angle) 0;这段代码中有几个技术细节值得注意使用real()函数处理复数结果避免绘图错误临界角后手动将透射率设为零虽然数学上反射率可能略高于1由于复数运算但物理上应该限制R≤1全反射现象在光纤通信中有重要应用正是利用这一原理光才能在光纤中长距离传输而不泄露能量。4. 可视化技巧与结果分析为了让计算结果更直观我通常采用多子图对比展示figure(Position,[100,100,1000,600]) % 振幅对比 subplot(2,2,1) plot(theta_i,abs(rs),b-,theta_i,abs(rp),r--,LineWidth,2) title(振幅反射系数) legend(s偏振,p偏振) xlabel(入射角(度)); ylabel(|r|) grid on % 能流对比 subplot(2,2,2) plot(theta_i,Rs,b-,theta_i,Rp,r--,theta_i,RsRp,k:,LineWidth,2) title(能流反射率) legend(R_s,R_p,R_sR_p) xlabel(入射角(度)); ylabel(反射率) grid on % 透射对比 subplot(2,2,3) plot(theta_i,Ts,b-,theta_i,Tp,r--,theta_i,TsTp,k:,LineWidth,2) title(能流透射率) legend(T_s,T_p,T_sT_p) xlabel(入射角(度)); ylabel(透射率) grid on % 能量守恒验证 subplot(2,2,4) plot(theta_i,RsTs,b-,theta_i,RpTp,r--,LineWidth,2) title(能量守恒验证) legend(s偏振,p偏振) xlabel(入射角(度)); ylabel(RT) ylim([0.99 1.01]) % 放大查看微小误差 grid on从这些曲线中可以观察到几个关键现象布儒斯特角处p偏振反射率为零临界角后全反射发生透射率为零能量守恒验证曲线应严格为1数值计算中的微小误差通常来自浮点运算精度在实际项目中我发现使用MATLAB的符号计算工具箱可以进一步提高计算精度特别是处理复数运算时% 使用符号计算提高精度 syms theta_i_sym r_p_sym (n2*cos(theta_i_sym) - n1*sqrt(1-(n1/n2*sin(theta_i_sym))^2)) / ... (n2*cos(theta_i_sym) n1*sqrt(1-(n1/n2*sin(theta_i_sym))^2)); % 可以精确计算布儒斯特角 brewster_angle atand(n2/n1);通过这些可视化分析我们不仅验证了菲涅尔公式的正确性还能直观理解不同偏振光在界面处的行为差异。这种理解对于设计光学系统、分析实验现象都具有重要价值。