jiexifa.m

本程序是用matlab仿真软件和解析法来计算光子晶体带隙的

% 解析法求解光子晶体带隙
% e1,e2介电常数，n1,n2折射率，c光速，d晶格常数，a,b薄膜厚度，u占空比，x相对频率，w0最大相对频率，point取点数，T解析矩阵；
clear;
e1=13;e2=1;u=0.5;w0=1.5;
n1=sqrt(e1);n2=sqrt(e2);
f=n1/n2;
d=0.1;a=u*d;b=d-a;
point=1000;
% x=w/c=2*pi*wi/d;
% wi=d/lambda;
% Ax=exp(i*x*n1*a);
% Bx=exp(i*x*n2*b);
% T=[1 1 -1 -1;-1 1 f -f;y/Ax y*Ax -Bx -1/Bx;-y/Ax y*Ax f*Bx -f/Bx];
% det(T)==0;
% 解方程 (4*Ax*f*Bx+Ax^2*f^2*y-Ax^2*f^2*y*Bx^2-2*Ax^2*f*y-2*Ax^2*f*y*Bx^2+y*Ax^2-y*Ax^2*Bx^2-y*f^2+y*f^2*Bx^2-2*y*f-2*y*f*Bx^2+4*y^2*f*Ax*Bx-y+y*Bx^2)/Ax/Bx=0;
% y=exp(i*k*d);
% y1=solve('(4*Ax*f*Bx+Ax^2*f^2*y-Ax^2*f^2*y*Bx^2-2*Ax^2*f*y-2*Ax^2*f*y*Bx^2+y*Ax^2-y*Ax^2*Bx^2-y*f^2+y*f^2*Bx^2-2*y*f-2*y*f*Bx^2+4*y^2*f*Ax*Bx-y+y*Bx^2)/Ax/Bx=0');
dx=w0*2*pi/d/point;
x=dx;
M=zeros(3,point);
for n=1:point
    Ax=exp(i*x*n1*a);
    Bx=exp(i*x*n2*b);
    t0=1/8/Ax/f/Bx;
    t1=-Ax^2*f^2+Ax^2*f^2*Bx^2+2*Ax^2*f+2*Ax^2*f*Bx^2-Ax^2+Ax^2*Bx^2+f^2-f^2*Bx^2+2*f+2*f*Bx^2+1-Bx^2;
    t2=(1+4*f+4*Ax^2*f^2+4*Ax^2*Bx^2+4*f^2*Bx^2-2*Ax^2+6*f^2-2*Bx^2-40*Ax^2*f^2*Bx^2-2*Ax^2*f^4-4*Ax^4*f+Ax^4*f^4-4*Ax^4*f^3+6*Ax^4*f^2+f^4*Bx^4-4*f^3*Bx^4+6*f^2*Bx^4-2*Ax^2*Bx^4-2*f^4*Bx^2+4*Ax^4*f*Bx^4+Ax^4+f^4+4*f^3+Bx^4-2*Ax^4*Bx^2+Ax^4*Bx^4+4*Ax^4*f^3*Bx^4+6*Ax^4*f^2*Bx^4-2*Ax^2*f^4*Bx^4+4*Ax^2*f^2*Bx^4-2*Ax^4*f^4*Bx^2+4*Ax^4*f^2*Bx^2+4*Ax^2*f^4*Bx^2+Ax^4*f^4*Bx^4-4*f*Bx^4)^(1/2);
    t3=t0*t1;
    M(2,n)=t3+t0*t2;
    M(3,n)=2*t3-M(2,n);
    M(1,n)=x;
    x=x+dx;
end
M
m=1;
for n=1:point
    if abs(real(M(2,m)))>1
       M(:,m)=[] ;
       m=m;
   else 
       m=m+1;
   end
end
M1=[M(1,:) M(1,:);M(2,:) M(3,:)];
n=1:2*(m-1);
figure
plot(real(-i*log(M1(2,n))/d),M1(1,n)*d/2/pi,'.')
axis([-pi/d,pi/d,0,w0])
title(['Photonic Bandgap  (n1=',num2str(n1),',n2=',num2str(n2),',d=',num2str(d),',u=',num2str(u),',w0=',num2str(w0),')'])
xlabel('k(-\pi/d<k<\pi/d)')
ylabel('d/\lambda(=\omegad/2\pic)')