MATLAB多层反射计算代码绘图仅输出单一值问题求助

问题原因&修正方案 #

核心故障点如下：

• 最直接的原因是代码中y_1变量未定义，Final= y_1*M_TOT;这行用到的y_1没有做任何初始化赋值，如果你的工作区之前遗留了固定值的y_1，会导致所有波长的计算结果完全相同，最终输出均匀直线。
• 其次反射系数的提取逻辑不符合传输矩阵的计算规则，不需要额外引入y_1做乘法，可以直接从总传输矩阵的元素提取反射系数。
• 层厚参数存在冗余：你的d数组第一个元素是2mm入射侧空气厚度，入射介质不需要计算传输相位，计入计算会引入不必要的误差。

修正后可运行代码 #

%This program is for calculating the total reflection from multilayer of
%epidermis-dermis-blood by considering air as the initial medium
clc;clear; % 加上clear清除工作区残留变量，避免之前的变量干扰
l=(0.3:0.03:3).*1.0E-3;%milimeter
e_inf=[1.0,3.0,5.5,2.1,1.0]; %Air-Epidermis-Dermis-Whole Blood_air(for reflection)
e_1=[1.0,58,58,130,1.0];
e_2=[1.0,3.6,5.6,3.8,1.0];
t_1=[0.1,9.4,9.4,14.4,0.1]*(1.0E-12);
t_2=[0.1,0.18,0.18,0.1,0.1]*(1.0E-12);
d=[0.2,0.002,0.05]*1.0E-3; % 去掉了冗余的入射空气层厚度，保留表皮、真皮、血液三层厚度
Q=zeros(length(l),5);

for x=1:1:length(l)
    M_TOT=eye(2,2);% 总矩阵初始化为单位矩阵，比ones更合理
    l_0=l(x); %wavelength
    w=((2*pi)*(3.0E8))./l_0; %angular frequency
    for z=2:1:4 % 从第二层（表皮）开始计算，对应介质层序号2-4，共3层
        Z=z+1; % 下一层介质序号
        % 计算当前层折射率
        epsilon_r_z=e_inf(1,z)+((e_1(1,z)-e_2(1,z))./(1+(w.*t_1(1,z)).^2))...
            +((e_2(1,z)-e_inf(1,z))./(1+(w.*t_2(1,z)).^2));
        epsilon_im_z=((e_1(1,z)-e_2(1,z)).*(w.*t_1(1,z)))./(1+(w.*t_1(1,z)).^2)...
            +((e_2(1,z)-e_inf(1,z)).*((w.*t_2(1,z)))./(1+(w.*t_2(1,z)).^2));
        n_1_z=sqrt((sqrt((epsilon_r_z).^2+(epsilon_im_z).^2)+epsilon_r_z)./2);
        n_2_z=sqrt((sqrt((epsilon_r_z).^2+(epsilon_im_z).^2)-epsilon_r_z)./2);
        n_z=n_1_z-1j*n_2_z; 

        % 计算下一层折射率
        epsilond_r_Z=e_inf(1,Z)+((e_1(1,Z)-e_2(1,Z))./(1+(w.*t_1(1,Z)).^2))...
            +((e_2(1,Z)-e_inf(1,Z))./(1+(w.*t_2(1,Z)).^2));
        epsilond_im_Z=((e_1(1,Z)-e_2(1,Z)).*(w.*t_1(1,Z)))./(1+(w.*t_1(1,Z)).^2)...
            +((e_2(1,Z)-e_inf(1,Z)).*((w.*t_2(1,Z)))./(1+(w.*t_2(1,Z)).^2));
        n_1_Z=sqrt((sqrt((epsilond_r_Z).^2+(epsilond_im_Z).^2)+epsilond_r_Z)./2);
        n_2_Z=sqrt((sqrt((epsilond_r_Z).^2+(epsilond_im_Z).^2)-epsilond_r_Z)./2);
        n_Z=n_1_Z-1j*n_2_Z; 

        phi_z=(n_z).*(2*pi/l_0).*d(1,z-1);% 取对应介质层的厚度
        r_TE_21_z=(n_Z-n_z)./(n_Z+n_z); 

        M_TE_z=[exp(-2j*phi_z) r_TE_21_z; r_TE_21_z.*exp(-2j*phi_z) 1]; 
        M_TOT=M_TE_z*M_TOT; % 保持左乘的规则不变
    end
    % 直接从总矩阵提取反射系数，删除y_1相关冗余代码
    r_12 = M_TOT(2,1)/M_TOT(1,1);
    Ref=abs(r_12).^2;     
    
    Q(x,1)=l_0*1.0E3;
    Q(x,2)=Ref;
end

plot(Q(:,1),Q(:,2))
title('Multilayer Reflection','Fontsize',14)
set(gcf,'color','w');
ylabel('Reflection','Fontsize',14)
xlabel('\lambda(mm)','Fontsize',14)

内容的提问来源于stack exchange，提问作者photonics_student