完整程序:

 

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

%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);