ttm.m

该程序用TMM法可计算一维介质的透射特性.

clear

na=3.4;
nb=1.0 ;
a=1;
b=4;
c=0.001*pi/(na*a+nb*b);  

for m=1:1000
    m
   for l=1:25
   %R(l)=0;
   %W(l)=0;
   R(l)=unifrnd(-0.5,0.5);
   W(l)=2*exp(-78.5*R(l)^2);
   end
for A=1:5100
 l=1;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x1=xb*xa;
 l=l+2;
   ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x2=xb*xa;
  l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x3=xb*xa;
   l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x4=xb*xa;   
     z1=x4*x3*x2*x1;
  l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x5=xb*xa;
 l=l+2;
   ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x6=xb*xa;
  l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x7=xb*xa;
   l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x8=xb*xa;   
      z2=x8*x7*x6*x5;
   l=l+2;
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x9=xb*xa;  
   l=l+2;   
     ba=A*c*na*a*(1+W(l)*R(l));
     bb=A*c*nb*b*(1+W(l+1)*R(l+1));
     d1=cos(ba);d2=sin(ba);
     d3=cos(bb);d4=sin(bb);
     xa=[d1,d2/na;-na*d2,d1];
     xb=[d3,d4/nb;-nb*d4,d3];
     x10=xb*xa;     
  z=x10*x9*z2*z1;
    %z=xa*z;
  % r=((z(2,2)-z(1,1))-i*(z(1,2)+z(2,1)))/((z(2,2)+z(1,1))-i*(z(1,2)-z(2,1)));
   %R(A)=abs(r)^2;
    t=2/((z(2,2)+z(1,1))-i*(z(1,2)-z(2,1)));
    T(A,m)=abs(t)^2;
end
end

for A=1:5100
    y(A)=0;
end
for A=1:5100
    for m=1:1000
        y(A)=y(A)+T(A,m);
    end
    Tp(A)=y(A)/1000;
end
s=linspace(1,5100,5100);
m=s*0.001;
q=Tp(s);
plot(m,q,'k')
axis([0  5  0  1.01])
xlabel('Ω');
ylabel('T');