clear all
wavelenth=600*10^(-9); %中心波長
N_w=1000;
n_a=4.6;n_b=1.6;n_base=1;   %介質(zhì)的折射率參數(shù)
d_a=0.25*wavelenth/n_a;  d_ai=(0.25-0.089)*wavelenth/n_a; d_a0=(0.25+0.089)*wavelenth/n_a;%A介質(zhì)的厚度
d_b=0.25*wavelenth/n_b;  d_bi=(0.25-0.089)*wavelenth/n_b; d_b0=(0.25+0.089)*wavelenth/n_b;%B介質(zhì)的厚度
for u=1:2          %TE和TM模的選擇
    for p=1:91         %角度值范圍
    theta(p)=0+(p-1).*1;    %入射角
    theta_0=theta(p)/180*pi;  %入射角弧度
        for t=1:(N_w+1)        % 頻率的變化范圍
            c=3.*10^(8);         %光速
            w_0=2.*pi.*c./wavelenth;  %中心波長的角頻率
            w_r=0+(t-1).*2./N_w;      %相對頻率的份數(shù)從0：2倍
            w=w_0.*w_r;                          %相對頻率的角頻率
            theta_a=asin(sin(theta_0)./n_a);    % A材料的折射角  
            theta_b=asin(sin(theta_0)./n_b);   %B材料折射角
            itate_a=n_a.*cos(theta_a);          % TE A導(dǎo)納
            dirta_a=w.*n_a.*d_a.*cos(theta_a)./c;dirta_ai=w.*n_a.*d_ai.*cos(theta_a)./c;dirta_a0=w.*n_a.*d_a0.*cos(theta_a)./c;%TE A位相角 
            itatm_a=n_a./cos(theta_a);          %TM  A導(dǎo)納
            itate_b=n_b.*cos(theta_b);           %TE  B導(dǎo)納
            dirta_b=w.*abs(n_b).*d_b.*cos(theta_b)./c;dirta_bi=w.*n_b.*d_bi.*cos(theta_b)./c;dirta_b0=w.*n_b.*d_b0.*cos(theta_b)./c;%TE  B位相角  
            itatm_b=n_b./cos(theta_b);            %TM B導(dǎo)納  
            ita0=[cos(theta_0),1/cos(theta_0)];%空氣介質(zhì)TE和TM模的選擇的導(dǎo)納計(jì)算方法
            ita_0=ita0(u);   %空氣介質(zhì)的導(dǎo)納
            theta_base=asin(sin(theta_0)./n_base);%入射角空氣
            itabase=[n_base.*cos(theta_base),n_base./cos(theta_base)];  %空氣介質(zhì)的導(dǎo)納
            ita_base=itabase(u);      %空氣介質(zhì)的U導(dǎo)納
            itaa=[itate_a,itatm_a];   %A介質(zhì)的U導(dǎo)納
            itab=[itate_b,itatm_b];   %B介質(zhì)的U導(dǎo)納             
            ita_a=itaa(u);ita_b=itab(u);%介質(zhì)的U導(dǎo)納選擇
            M_A=[cos(dirta_a),-i.*sin(dirta_a)./ita_a;-i.*sin(dirta_a).*ita_a,cos(dirta_a)];%A介質(zhì)的矩陣
            M_B=[cos(dirta_b),-i.*sin(dirta_b)./ita_b;-i.*sin(dirta_b).*ita_b,cos(dirta_b)];%B介質(zhì)的矩陣
            M_Ai=[cos(dirta_ai),-i.*sin(dirta_ai)./ita_a;-i.*sin(dirta_ai).*ita_a,cos(dirta_ai)];%Ai介質(zhì)的矩陣
            M_Bi=[cos(dirta_bi),-i.*sin(dirta_bi)./ita_b;-i.*sin(dirta_bi).*ita_b,cos(dirta_bi)];%Bi介質(zhì)的矩陣
            M_A0=[cos(dirta_a0),-i.*sin(dirta_a0)./ita_a;-i.*sin(dirta_a0).*ita_a,cos(dirta_a0)];%Ai介質(zhì)的矩陣
            M_B0=[cos(dirta_b0),-i.*sin(dirta_b0)./ita_b;-i.*sin(dirta_b0).*ita_b,cos(dirta_b0)];%Bi介質(zhì)的矩陣
            M_total=(M_Bi*M_Ai)^8*(M_B*M_A)^10*(M_A*M_B)^10*(M_A0*M_B0)^8;
            %M_total=(M_A0*M_B0)^8;
            A=M_total(1,1);B=M_total(1,2);C=M_total(2,1);D=M_total(2,2);                    %矩陣個元素的值
            R(t,p,u)=abs((A.*ita_0+B.*ita_0.*ita_base-C-D.*ita_base)./(A.*ita_0+B.*ita_0.*ita_base+C+D.*ita_base)).^2;%反射率
            T(t,p,u)=abs(2.*ita_0./(A.*ita_0+B.*ita_0.*ita_base+C+D.*ita_base)).^2;  %透射率
        end
    end
end
T_te=T(:,:,1);T_tm=T(:,:,2);T_total=[fliplr(T_tm(:,2:91)),T_te];
R_te=R(:,:,1);R_tm=R(:,:,2);R_total=[fliplr(R_tm(:,2:91)),R_te];
clear j;clear t;clear k; clear p; clear u;clear s;