bandgapte.m

calculating 1D photonic crystal

clear all;
omega=linspace(0,12*pi*10^14,1000);

per=10^-6;c=3*10^8;
% beta=linspace(10^-7,6*pi*10^-6,1000);
 beta=linspace(0,0,1000);
n1=1.5;
n2=2.2;
a=3*10^-7;
b=per-a;
omega0=omega.*per./c;
beta0=beta.*per;
k1x=sqrt((omega.*n1/c).^2-beta.^2);
k2x=sqrt((omega.*n2/c).^2-beta.^2);
A=exp(-i*k1x.*a).*(cos(k2x.*b)-0.5*i*(k2x./k1x+k1x./k2x).*sin(k2x.*b));
B=exp(i*k1x.*a).*(-0.5*i*(k2x./k1x-k1x./k2x).*sin(k2x.*b));
C=exp(-i*k1x.*a).*(0.5*i*(k2x./k1x-k1x./k2x).*sin(k2x.*b));
D=exp(i*k1x.*a).*(cos(k2x.*b)+0.5*i*(k2x./k1x+k1x./k2x).*sin(k2x.*b));
    K=acos((A+D)./2)./per;
    hold on;
subplot(2,1,1),plot(omega0,real(K),'b');

subplot(2,1,2),plot(omega0,abs(imag(K)),'r');
xte=zeros(1,1000);
zte=zeros(1,1000);
m=1;
for j=1:1:1000
    if abs(imag(K(j)))>0
    xte(1,m)=beta0(1,j);zte(1,m)=omega(1,j);
    m=m+1;
end;
end;
% % axis([0+m,2+m*2,m*pi*10^6,(m+1)*pi*10^6]);
% 
% % end;
% % plot(K,A+D)
% if abs(A+D)<2
% plot(beta0,omega0,'r');
% end;