Need Urgent Help:Matlab code for Transmission coefficient for a Pc

[elham.n]

Hi,
I have written MATLAB code for Transmission coeficient for a one-dimensional photonic crystal,but i have problems with it.it does not result the right answer!
could anyone help me about this?
i use Transfer matrix method:
An example:
an ENG-DPS structure:


%Transmission for a multilayer structure versus frequency(GHz)

clc; clear all; close all
c0=3*10^8;
d=0.01; %meter
dE=1.5*d;
dP=0.5*d;
N=20; %N is the number of periods.
w0=pi*c0/d;
wep=sqrt(2)*w0;
epsP=1; muP=1;
muE=1;
teta=0*pi/180;

%%%%%%%%%%%%%%%%%%%%%

for j=1:1501;
w(j)=0.5*w0+(j-1)*0.001*w0;
epsE(j)=1-(wep^2)/(w(j)^2);
qE(j)=(sqrt(epsE(j))/sqrt(muE))*sqrt(1-(sin(teta))^2/epsE(j)*muE);
qP=(sqrt(epsP))/sqrt(muP)*sqrt(1-(sin(teta))^2/epsP*muP);
kE(j)=(w(j)/c0)*sqrt(epsE(j))*sqrt(muE)*sqrt(1-(sin(teta))^2/epsE(j)*muE);
kP(j)=(w(j)/c0)*sqrt(epsP)*sqrt(muP)*sqrt(1-(sin(teta))^2/epsP*muP);
ME=[cos(kE(j)*dE) i*(1/qE(j))*sin(kE(j)*dE); i*qE(j)*sin(kE(j)*dE) cos(kE(j)*dE)];
MP=[cos(kP(j)*dP) i*(1/qP)*sin(kP(j)*dP); i*qP*sin(kP(j)*dP) cos(kP(j)*dP)];
M=(ME*MP)^N;
t(j)=2*cos(teta)/((M(2,2)*cos(teta)+M(1,1)*cos(teta))-(M(2,1)+M(1,2)*cos(teta)*cos(teta)));
T(j)=(abs(t(j)))^2;
end

%%%%%%%%%%%%%%%%

%figure

plot(w/w0,T,'b','linewidth',1.5)
set(gca,'fontsize',15);
xlabel('\omega/\omega_0');
ylabel('Transmission');
axis([0,2,0,1.1]);
title(['ENG-Air',', N=',num2str(N),', \theta=',num2str(teta*180/pi),],'FontSize',12 );
hold on

 

[NateD]

%%%%%%%%%%%%%%%for j=1:1501;w(j)=0.5*w0+(j-1)*0.001*w0;epsE(j)=1-(wep^2)/(w(j)^2);qE(j)=(sqrt(epsE(j))/sqrt(muE))*sqrt(-1-(sin(teta))^2/epsE(j)*muE);kE(j)=(w(j)/c0)*sqrt(epsE(j))*sqrt(muE)*sqrt(-1-(sin(teta))^2/epsE(j)*muE);ME=[cos(kE(j)*dE) i*(1/qE(j))*sin(kE(j)*dE); i*qE(j)*sin(kE(j)*dE) cos(kE(j)*dE)];M=(ME*MP)^N;t(j)=2*cos(teta)/((M(2,2)*cos(teta)+M(1,1)*cos(teta))-(M(2,1)+M(1,2)*cos(teta)*cos(teta)));T(j)=(abs(t(j)))^2;end%%%%%%%%%%%%%%%%%%%%%%%%plot(w/w0,T,'r','linewidth',1.5)hold off%%%%%%%%%%%%%%%%%%%%%%%The result is wrong. the transmission coefficient should be approximately zero at the bandgap frequencies. but in my case, the transmission coefficient has a peak near the bandgap frequency!thanks in advance. </code>A:There is a mistake in the calculation of <code>qE</code> in the second loop:<code>qE(j)=(sqrt(epsE(j))/sqrt(muE))*sqrt(-1-(sin(teta))^2/epsE(j)*muE);</code>The minus sign in front of <code>(sin(