基于matlab实现的光折射反射（不同界面）程序

完整程序:

%平面电磁波在不同介质界面上入射、反射、折射仿真

%ReadMe!!!在下述说明的用户输入区内输入入射角和两介质折射率，
%输出反射折射示意图与反射折射系数随入射角变化的曲线
%—————————————————————————————————————

%——————————用户输入区起始——————————%
theta1 = 30; %请输入入射角（大于0且小于等于90°）
n1 = 1.5; %请输入介质1的折射率
n2 = 1; %请输入介质2的折射率
%——————————用户输入区结束——————————%

flag = 0; %表示仿真状态，0代表正常入射，1代表垂直入射，2代表全反射
format short
if theta1 == 0 %垂直入射
flag = 1;
elseif sind(theta1)*n1/n2 >= 1 %满足全反射条件
flag = 2;
else %其他通常情况的折射反射
flag = 0;
end

if flag == 0 %当不满足全反射条件时
theta2 = asind(sind(theta1)*n1/n2); %求出折射角theta2
rs = -sind(theta1-theta2)/sind(theta1 + theta2); %s波反射系数
ts = 2*sind(theta2)*cosd(theta1)/sind(theta1+theta2); %s波折射系数
rp = tand(theta1-theta2)/tand(theta1+theta2); %p波反射系数
tp = 2*sind(theta2)*cosd(theta1)/(sind(theta1+theta2)*cosd(theta1-theta2)); %p波折射系数
elseif flag == 2 %当满足全反射条件时
rs=1;rp=1;tp=3;ts=2;
elseif flag == 1 %当垂直入射
theta2 = 0;
rs=-1;rp=-1;tp=0;ts=0;
end

figure('NumberTitle', 'off', 'Name', '图一反射折射示意图'); %绘制光在两介质表面的反射、折射示意图
r = 5; %绘制时的线长

plot([-3 3],[0 0],'k--','LineWidth',1); %绘制辅助线
hold on;
plot([0 0],[-2 2],'k--','LineWidth',1);
hold on;

if flag ~= 1 %判断是否为垂直入射
x = -sind(theta1)*r:0.1:0;
y = -x./tand(theta1);
elseif flag == 1
y = 0:0.1:r;
x = -0.03*ones(size(y));
end
RuShe = plot(x,y,'r','LineWidth',2); %绘制入射光线
hold on;
FanShe = plot(-x,y,'g','LineWidth',2); %绘制反射光线
legend([RuShe FanShe],'入射光','反射光');
x_scale = max(abs(r*sind(theta1)), 4.5); %确定坐标轴范围
y_scale = max(abs(r*cosd(theta1)), 4.5);
if flag ~= 2 %判断是否为全反射
hold on;
if flag ~= 1 %判断是否为垂直入射
x = 0:0.1:r*sind(theta2);
y = -x./tand(theta2);
else
y = -y;
end
ZheShe = plot(x,y,'b','LineWidth',2);
legend([RuShe FanShe ZheShe],'入射光','反射光','折射光');
x_scale = max(x_scale,abs(r*sind(theta2))); %加入折射光线后，重新确定坐标轴范围
y_scale = max(y_scale,abs(r*cosd(theta2)));
end

grid on;
axis equal;
axis ([-x_scale x_scale -y_scale y_scale]);
title('折射，反射示意图','FontSize',12);

%在图中添加说明
if flag ~= 2
temp = num2str(theta2);
else
temp = '全反射';
end
textstr = {['n1=' num2str(n1) ', n2=' num2str(n2)]; ['入射角为： ' num2str(theta1) '°']; ['折射角为： ' temp '°']; ['rs = ' num2str(rs)]; ['rp = ' num2str(rp)]; ['ts = ' num2str(ts)]; ['tp = ' num2str(tp)]};
text(min(-x_scale+0.5, -4), -1.8, textstr, 'FontSize', 10);

figure('NumberTitle', 'off', 'Name', '图二反射折射系数随入射角的变化曲线'); %绘制各反射、折射系数随入射角变化的曲线图
theta_1 = linspace(0, 90, 300); %入射角
theta_2 = asind(sind(theta_1)*n1/n2); %求出折射角theta2
rs = -sind(theta_1-theta_2)./sind(theta_1 + theta_2); %s波反射系数
ts = 2*sind(theta_2).*cosd(theta_1)./sind(theta_1+theta_2); %s波折射系数
rp = tand(theta_1-theta_2)./tand(theta_1+theta_2); %p波反射系数
tp = 2*sind(theta_2).*cosd(theta_1)./(sind(theta_1+theta_2).*cosd(theta_1-theta_2)); %p波折射系数
if n1 > n2 %当可能发生全反射时
theta_c = asind(n2/n1); %全反射角
temp = find(theta_1 >= theta_c);
rs(temp)=1; rp(temp)=1; tp(temp)=3; ts(temp)=2;
end
Rs = plot(theta_1, rs, 'r', 'LineWidth', 1.5);
hold on
Ts = plot(theta_1, ts, 'g', 'LineWidth', 1.5);
hold on
Rp = plot(theta_1, rp, 'b', 'LineWidth', 1.5);
hold on
Tp = plot(theta_1, tp, 'y', 'LineWidth', 1.5);
hold on
plot([0 90], [0 0], 'k', 'LineWidth', 0.5); %辅助线
axis([0 90 -1 1]);
if n1 > n2
hold on
plot([0 theta_c], [1 1], 'k--', 'LineWidth', 0.5); %辅助线
hold on
plot([0 theta_c], [2 2], 'k--', 'LineWidth', 0.5);
hold on
plot([theta_c theta_c], [-1 3.5], 'k--', 'LineWidth', 0.5);
axis([0 90 -1 3.5]);
end
grid on;
legend([Tp Ts Rs Rp], 't_{p}', 't_{s}', 'r_{s}', 'r_{p}');
title('t_{p},t_{s},r_{s},r_{p}随\theta_{1}的变化', 'FontSize', 12);

figure('NumberTitle', 'off', 'Name', '图三反射比透射比曲线图');
rhos=(sind(theta_1-theta_2)).^2./(sind(theta_1+theta_2)).^2;
taus=4*n2*cosd(theta_2).*(sind(theta_2)).^2.*(cosd(theta_1)).^2./(n1*cosd(theta_1).*sind(theta_1+theta_2).^2);
rhop=(tand(theta_1-theta_2)).^2./(tand(theta_1+theta_2)).^2;
taup=4*n2*cosd(theta_2).*(sind(theta_2)).^2.*(cosd(theta_1)).^2./(n1*cosd(theta_1).*sind(theta_1+theta_2).^2.*cosd(theta_1-theta_2).^2);

if n1 > n2 %当可能发生全反射时
rhos(temp)=1; rhop(temp)=1; taup(temp)=0; taus(temp)=0;
end
Rhos = plot(theta_1, rhos, 'r', 'LineWidth', 1.5);
hold on
Taus = plot(theta_1, taus, 'g', 'LineWidth', 1.5);
hold on
Rhop = plot(theta_1, rhop, 'b', 'LineWidth', 1.5);
hold on
Taup = plot(theta_1, taup, 'y', 'LineWidth', 1.5);
hold on
if n1 > n2
hold on
plot([theta_c theta_c], [0 1], 'k--', 'LineWidth', 0.5);
end
axis([0 90 0 1.5]);
grid on;
legend([Taup Taus Rhos Rhop], '{\tau}_{p}', '{\tau}_{s}', '{\rho}_{s}', '{\rho}_{p}');
title('\tau_{p},\tau_{s},\rho_{s},\rho_{p}随\theta_{1}的变化', 'FontSize', 12);