fzs.m

matlat程序用于计算带金属覆层的一维光子晶体的反射率

clc,clear
a=59.0;                   
b=115.0;    
hAg=42;
special=39.0;

nb=1.45;          
hba=6.626196/10^34/2/pi;
N=9;          
c=299792458;           
theta=20/180*pi;
lambda=[206.6 213.8 221.4 229.6 238.4 248.0 253.0 258.3 263.8 269.5 275.5 281.8 288.3 295.2 298.8 302.4 306.1 310.0 311.5 313.9 315.5 317.9 319.5 322.0 323.7 326.3 330.6 332.4 335.1 339.7 344.4 354.2 364.7 375.7 387.5 400.0 413.3 427.5 442.8 459.2 476.9 495.9 516.6 539.1 563.6 590.4 619.6 652.6 688.8 729.3 774.9 826.6 885.6 953.7 1033 1127 1240];
n1=[1.125 1.173 1.208 1.238 1.265 1.298 1.320 1.343 1.372 1.404 1.441 1.476 1.502 1.519 1.522 1.496 1.432 1.323 1.246 1.149 1.044 0.932 0.815 0.708 0.616 0.526 0.371 0.321 0.294 0.259 0.238 0.209 0.186 0.200 0.192 0.173 0.173 0.160 0.157 0.144 0.132 0.130 0.130 0.129 0.120 0.121 0.131 0.140 0.140 0.148 0.143 0.145 0.163 0.198 0.226 0.251 0.329];                %57
n2=[1.27 1.29 1.30 1.31 1.33 1.35 1.35 1.35 1.35 1.33 1.31 1.26 1.19 1.08 0.992 0.882 0.766 0.647 0.586 0.540 0.514 0.504 0.526 0.565 0.609 0.663 0.813 0.902 0.986 1.12 1.24 1.44 1.61 1.67 1.81 1.95 2.1 1 2.26 2.40 2.56 2.72 2.88 3.07 3.25 3.45 3.66 3.88 4.15 4.44 4.74 5.09 5.50 5.95 6.43 6.99 7.67 8.49];
p=450.1:.1:850;
for t=1:4000   
  for z=1:1:55
      j=z+1;
    if lambda(z)<=p(t)&lambda(j)>=p(t)
       nAg(t)=n1(z)+(p(t)-lambda(z))/(lambda(j)-lambda(z))*(n1(j)-n1(z))+i*(n2(z)+(p(t)-lambda(z))/(lambda(j)-lambda(z))*(n2(j)-n2(z)));
    end
  end
  omega=2*pi*c/p(t);
Xi(t)=1-6.2^2/(hba^2*(omega*10^9/1.60217646*10^19)^2-3.8^2+i*5.0*10^(-4)*hba*omega*10^9/1.60217646*10^19);
na=Xi(t)^0.5;
  

  x=omega*nAg(t)*cos(theta)/c;;
  M3=[cos(x*hAg),-i*sin(x*hAg)/cos(theta)/nAg(t);-i*cos(theta)*nAg(t)*sin(x*hAg),cos(x*hAg)];
   alpha=omega*na*cos(theta)/c;
   beta =omega*nb*cos(theta)/c;
  Ma=[cos(alpha*a),-i*sin(alpha*a)/cos(theta)/na;-i*cos(theta)*na*sin(alpha*a),cos(alpha*a)];                   %求出频率为w时介质层na的传输矩阵
  Mb=[cos(beta*b),-i*sin(beta*b)/cos(theta)/nb;-i*cos(theta)*nb*sin(beta*b),cos(beta*b)];                      %求出频率为w时介质层nb的传输矩阵
  M1=[cos(alpha*special),-i*sin(alpha*special)/cos(theta)/na;-i*cos(theta)*na*sin(alpha*special),cos(alpha*special)]*Mb;

  Mnon=M1*((Ma*Mb)^N);
  rnon=((Mnon(1,1)*cos(theta)+Mnon(1,2)*cos(theta)*cos(theta)-Mnon(2,1)-Mnon(2,2)*cos(theta))/(Mnon(1,1)*cos(theta)+Mnon(1,2)*cos(theta)*cos(theta)+Mnon(2,1)+Mnon(2,2)*cos(theta)));
  Rnon(t)=abs(rnon)^2;
  M=M3*M1*((Ma*Mb)^N);                                          %求出总的传输矩阵
  r=((M(1,1)*cos(theta)+M(1,2)*cos(theta)*cos(theta)-M(2,1)-M(2,2)*cos(theta))/(M(1,1)*cos(theta)+M(1,2)*cos(theta)*cos(theta)+M(2,1)+M(2,2)*cos(theta)));   
  R(t)=abs(r)^2;%求出反射率
end 
plot(p,R,p,Rnon,':'),xlabel('\lambda'),ylabel('R');


%subplot(2,2,1)
%plot(p,R)
%xlim([685,779])
%subplot(2,2,2)
%plot(p,Rnon)
%subplot(2,2,3)
%plot(p,abs(nAg))
%xlim([685,779])