MATLAB利用菲涅尔公式仿真光的折射

从空气进入石英玻璃，计算振幅反射率、透射率和对应的绝对值：

clear;close all;n1 = 1; %空气折射率n2 = 1.45;%平板玻璃折射率theta = 0:0.1:90;%角度从0到90a = theta*pi/180;%角度转化为弧度rp = (n2*cos(a)-n1*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));tp = (2*n1*cos(a))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));rs = (n1*cos(a)-n2*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));ts = (2*n1*cos(a))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));figure(1);subplot(1,2,1);plot(theta,rp,'-',theta,rs,'--',theta,abs(rp),':' ,theta,abs(rs),'-.','linewidth',2);legend('r_p','r_s','|r_p|','|r_s|');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 -1 1]);grid on;subplot(1,2,2);plot(theta,tp,'-',theta,ts,'--',theta,abs(tp),':' ,theta,abs(ts),'-.','linewidth',2);legend('t_p','t_s','|t_p|','|t_s|');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 0 1]);grid on;

结果如下：

可以看出：

1、首先，当入射角等于0时，有反射波也有透射波，这个时候并不是完全进入折射率高的介质内；

2、当入射角等于90时，透射波为0，能量完全反射。

3、当角度增大时，投射的两个偏振都减小，反射的s偏振逐渐增大，p偏振先减小再增大。

4、t始终大于0，折射光的相位和入射光始终是一致的；

5、rs始终小于0，表明s偏振光在反射面上发生了π相位变化；

6、rp前期相位变化为0，后期相位相差π。

反之，从光密介质进入光疏介质，全反射存在

%从光密到光疏clear;close all;n1 = 1.45; %空气折射率n2 = 1;%平板玻璃折射率theta = 0:0.1:90;%角度从0到90a = theta*pi/180;%角度转化为弧度rp = (n2*cos(a)-n1*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));tp = (2*n1*cos(a))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));rs = (n1*cos(a)-n2*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));ts = (2*n1*cos(a))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));figure(1);subplot(1,2,1);plot(theta,rp,'-',theta,rs,'--',theta,abs(rp),':' ,theta,abs(rs),'-.','linewidth',2);legend('r_p','r_s','|r_p|','|r_s|');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 -1 1.2]);grid on;subplot(1,2,2);plot(theta,tp,'-',theta,ts,'--',theta,abs(tp),':' ,theta,abs(ts),'-.','linewidth',2);legend('t_p','t_s','|t_p|','|t_s|');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 0 3]);grid on;可以看出：1、对于rp仍然存在布鲁斯特角；2、光密介质进入光疏介质存在全反射现象，rp和rs的绝对值为0；3、之所以出现如此大的差异，是因为当角度超过全反射角时，rs、rp等为复数。matlab中plot函数在绘图时会忽视复数的虚部。【例1.4】对于光密到光疏介质，我们可以计算能流密度来表征：%从光密到光疏clear;close all;n1 = 1.45; %空气折射率n2 = 1;%平板玻璃折射率theta = 0:0.1:90;%角度从0到90a = theta*pi/180;%角度转化为弧度rp = (n2*cos(a)-n1*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));tp = (2*n1*cos(a))./(n2*cos(a)+n1*sqrt(1-(n1/n2)^2*(sin(a)).^2));rs = (n1*cos(a)-n2*sqrt(1-(n1/n2)^2*(sin(a)).^2))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));ts = (2*n1*cos(a))./(n1*cos(a)+n2*sqrt(1-(n1/n2)^2*(sin(a)).^2));Rp = abs(rp).^2;Rs = abs(rs).^2;Rn = (Rp+Rs)/2;Tp = n2*sqrt(1-(n1/n2)^2*(sin(a)).^2)./(n1*cos(a)).*abs(tp).^2;Ts = n2*sqrt(1-(n1/n2)^2*(sin(a)).^2)./(n1*cos(a)).*abs(ts).^2;Tn = (Tp+Ts)/2;subplot(1,2,1);plot(theta,Rp,'-',theta,Rs,'-.',theta,Rn,'--','linewidth',2);legend('R_p','R_s','R_n');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 0 1]);grid on;subplot(1,2,2);plot(theta,Tp,'-',theta,Ts,'-.',theta,Tn,'--','linewidth',2);legend('T_p','T_s','T_n');xlabel('\theta_i');ylabel('Amplitude');title(['n_1=',num2str(n1),',n_2=',num2str(n2)]);axis([0 90 0 1]);grid on;

结果如下：

可以看出：

1、对于rp仍然存在布鲁斯特角；

2、光密介质进入光疏介质存在全反射现象，rp和rs的绝对值为0；

3、之所以出现如此大的差异，是因为当角度超过全反射角时，rs、rp等为复数。matlab中plot函数在绘图时会忽视复数的虚部。

来源《高等光学仿真——光波导、激光》