bandgaptm.m

calculating 1D photonic crystal

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

per=10^-6;c=3*10^8;
namta=2*pi*c./omega;
% beta=linspace(10^-7,6*pi*10^-6,1000);
%  beta=linspace(0,0,1000);
theta=0;
beta=2*pi./namta.*sin(theta);
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*(n1^2.*k2x./k1x./n2^2+n2^2.*k1x./k2x./n1^2).*sin(k2x.*b));
B=exp(i*k1x.*a).*(0.5*i*(n1^2.*k2x./k1x./n2^2-n2^2.*k1x./k2x./n1^2).*sin(k2x.*b));
C=exp(-i*k1x.*a).*(-0.5*i*(n1^2.*k2x./k1x./n2^2-n2^2.*k1x./k2x./n1^2).*sin(k2x.*b));
D=exp(i*k1x.*a).*(cos(k2x.*b)+0.5*i*(n1^2.*k2x./k1x./n2^2+n2^2.*k1x./k2x./n1^2).*sin(k2x.*b));
%   for m=0:5
hold on;
    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')
% % 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;