分层介质反射和透射系数传播矩阵法

写论文用到了，找了好久，在电磁波理论（葛德彪，魏兵著）找到了相应的代码，发到网上吧，希望大家方便点。
第一个：分层介质反射和透射系数传播矩阵法

clear
epsle_0=8.854187817e-12;  %air epsion
mu0=4*pi*1e-7;                  %air mu0
epsion1=1;
epsion2=2;
epsion3=3;
% % epsion=epsle_0*[1,4.2,epsion1,4.2, epsion2,4.2,epsion3,4.2,1];%9
% % mur=mu0*[1,1,1,1,1,1,1,1,1];
% theta=90*pi/180;
%  sigma=[0,0,0,0,0,0,0,0,0];
%  sigmam=[0,0,0,0,0,0,0,0,0];
% d=[0,-0.001,-0.001006,-0.001012,-0.002012,-0.003012,-0.003018,-0.004018,-0.005018];


epsion=epsle_0*[1,4,1,4,1,4,1,4,1,4,1,4,1];%9
mur=mu0*[1,1.,1,1.,1,1.,1,1.,1,1.,1,1.,1];

d=[0,-0.01,-0.02,-0.03,-0.04,-0.05,-0.06,-0.07,-0.08,-0.09,-0.1,-0.11,-0.12];

 sigma=[0,0,0,0,0,0,0,0,0,0,0,0,0];
 sigmam=[0,0,0,0,0,0,0,0,0,0,0,0,0];

theta=30*pi/180;

M=size(epsion);
N=M(1,2)-1;

disp('Please choose the incident wave:TM=-1 or TE=1');
wavetype=input('wavetype is:');

fid1=fopen('RT.dat','wt');
for f=10:10:10000
     

   omega=2*pi*f*1e6;%入射波圆频率（实际频率为循环变量的1e6倍）
  k = omega*sqrt(epsion .* mur); 
     kx=k(:,1)*sin(theta);%入射波矢量在X方向的投影（连续）
     
     kz=sqrt(k.^2-kx.^2);%存放各介质区域的波矢量z分量大小矩阵
     
     
     if wavetype==1
         for n=1:N
            p(n)=mur(:,n+1)*kz(:,n)/(mur(:,n)*kz(:,n+1));
         end
     elseif wavetype==-1
         for n=1:N
             p(n)=epsion(:,n+1)*kz(:,n)/(epsion(:,n)*kz(:,n+1));
         end
     end
     
     for n=1:N
         r(n)=(1-p(:,n))/(1+p(:,n));
     end
     VV=eye(2);
     for n=1:N-1
         v=0.5*(1+p(:,n));
         V11=v*exp(-j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V12=v*r(:,n)*exp(-j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V21=v*r(:,n)*exp(j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V22=v*exp(j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V=[V11,V12;V21,V22];
         VV=VV*V;
     end
     
     v1=0.5*(1+p(:,N));
     Vtn11=v1*exp(j*kz(:,N+1)*d(:,N));
     Vtn12=v1*r(:,N)*exp(j*kz(:,N+1)*d(:,N));
     Vtn21=v1*r(:,N)*exp(-j*kz(:,N+1)*d(:,N));
     Vtn22=v1*exp(-j*kz(:,N+1)*d(:,N));
     Vtn=[Vtn11,Vtn12;Vtn21,Vtn22];
     
     R=(-VV(1,2)/VV(1,1))*exp(j*2*kz(:,1)*d(:,1));
     
     T=(-VV(2,1)*VV(1,2)/VV(1,1)+VV(2,2))*exp(-j*(kz(:,N+1)*d(:,N)-kz(:,1)*d(:,1)));
     
     plot(f,abs(R),'*',f,abs(T),'o');
     hold on
     
     fprintf(fid1,'\n%8d%+12.5E%+12.5E',f,abs(R),abs(T));
end

legend('反射系数','透射系数');
xlabel('频率 GHZ');
ylabel('反射透射系数');
fclose(fid1);    

TE波：
TM波：

下面是课本给出的

还有个分层介质反射系数的连分法
下面展示一些 内联代码片。

// A code block
var foo = 'bar';

%%%%%%分层反射连数法

clear
epsle_0=8.854187817e-12;  %air epsion
mu0=4*pi*1e-7;                  %air mu0
epsion1=1;
epsion2=2;
epsion3=3;
% % epsion=epsle_0*[1,4.2,epsion1,4.2, epsion2,4.2,epsion3,4.2,1];%9
% % mur=mu0*[1,1,1,1,1,1,1,1,1];
% theta=90*pi/180;
%  sigma=[0,0,0,0,0,0,0,0,0];
%  sigmam=[0,0,0,0,0,0,0,0,0];
% d=[0,-0.001,-0.001006,-0.001012,-0.002012,-0.003012,-0.003018,-0.004018,-0.005018];


epsion=epsle_0*[1,4,1,4,1,4,1,4,1,4,1,4,1];%9
mur=mu0*[1,1.,1,1.,1,1.,1,1.,1,1.,1,1.,1];

d=[0,-0.01,-0.02,-0.03,-0.04,-0.05,-0.06,-0.07,-0.08,-0.09,-0.1,-0.11,-0.12];

 sigma=[0,0,0,0,0,0,0,0,0,0,0,0,0];
 sigmam=[0,0,0,0,0,0,0,0,0,0,0,0,0];

theta=30*pi/180;



disp('Please choose the incident wave:TM=-1 or TE=1');
wavetype=input('wavetype is:');

fid1=fopen('RT.dat','wt');
for f=100:50:10000
     

   omega=2*pi*f*1e6;%入射波圆频率（实际频率为循环变量的1e6倍）
   
      k = omega*sqrt(epsion .* mur); 
     kx=k(:,1)*sin(theta);%入射波矢量在X方向的投影（连续）
     
     kz=sqrt(k.^2-kx.^2);%存放各介质区域的波矢量z分量大小矩阵
     M=size(epsion);
N=M(1,2)-1;
     
     if wavetype==1
         for n=1:N
            p(n)=mur(:,n+1)*kz(:,n)/(mur(:,n)*kz(:,n+1));
         end
     elseif wavetype==-1
         for n=1:N
             p(n)=epsion(:,n+1)*kz(:,n)/(epsion(:,n)*kz(:,n+1));
         end
     end
     
     for n=1:N
         r(n)=(1-p(:,n))/(1+p(:,n));
     end
     R=0;
         for n=N:-1:1
             A=1/r(:,n);
             B=exp(i*2*kz(:,n+1)*d(:,n));
             C=exp(i*2*kz(:,n)*d(:,n));
             R=A*C+(1-A^2)*B*C/(A*B+R);
         end
         plot(f,abs(R),'*');
         hold on;
     
     VV=eye(2);
     for n=1:N-1
         v=0.5*(1+p(:,n));
         V11=v*exp(-j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V12=v*r(:,n)*exp(-j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V21=v*r(:,n)*exp(j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V22=v*exp(j*kz(:,n+1)*(d(:,n+1)-d(:,n)));
         V=[V11,V12;V21,V22];
         VV=VV*V;
     end
end
%      v1=0.5*(1+p(:,N));
%      Vtn11=v1*exp(j*kz(:,N+1)*d(:,N));
%      Vtn12=v1*r(:,N)*exp(j*kz(:,N+1)*d(:,N));
%      Vtn21=v1*r(:,N)*exp(-j*kz(:,N+1)*d(:,N));
%      Vtn22=v1*exp(-j*kz(:,N+1)*d(:,N));
%      Vtn=[Vtn11,Vtn12;Vtn21,Vtn22];
%      
%      R=(-VV(1,2)/VV(1,1))*exp(j*2*kz(:,1)*d(:,1));
%      
%      T=(-VV(2,1)*VV(1,2)/VV(1,1)+VV(2,2))*exp(-j*(kz(:,N+1)*d(:,N)-kz(:,1)*d(:,1)));
%      
%      plot(f,abs(R),'*',f,abs(T),'o');
%      hold on
%      
%      fprintf(fid1,'\n%8d%+12.5E%+12.5E',f,abs(R),abs(T));
% end
% 
% legend('反射系数','透射系数');
% xlabel('频率 GHZ');
% ylabel('反射透射系数');
% fclose(fid1);
     

TE波
TM波：
课本给的：
课本链接：

通过百度网盘分享的文件：电磁波理论_葛德彪，魏兵著_978-7-03-032006-3.pdf
链接：https://pan.baidu.com/s/1YmZ3uOLFt9j9WiMW2a2_Dg
提取码：g8js