pwe1.m

2D photonic crystal Band Structure

function PWE(N,RES,r,epa,epb)
clc;
L = 1;
Band1 = zeros(1,3*RES+1);
Band2 = zeros(1,3*RES+1);
Band3 = zeros(1,3*RES+1);
Band4 = zeros(1,3*RES+1);
Band5 = zeros(1,3*RES+1);

Band1TM = zeros(1,3*RES+1);
Band2TM = zeros(1,3*RES+1);
Band3TM = zeros(1,3*RES+1);
Band4TM = zeros(1,3*RES+1);
Band5TM = zeros(1,3*RES+1);

index = 0;
% GAMMA-X section
for KappaX=0:pi/L/RES:pi/L
    KappaX
    index = index + 1;
    KappaY = 0;
    W   = PWETE(KappaX,KappaY,N,L,r,epa,epb);
    WTM = PWETM(KappaX,KappaY,N,L,r,epa,epb);
    
    dTM = eig(WTM);
    dTM = sort(sqrt(dTM));
    Band1TM(index) = dTM(1);
    Band2TM(index) = dTM(2);
    Band3TM(index) = dTM(3);
    Band4TM(index) = dTM(4);
    Band5TM(index) = dTM(5);
    
    d = eig(W);
    d = sort(sqrt(d));
    Band1(index) = d(1);
    Band2(index) = d(2);
    Band3(index) = d(3);
    Band4(index) = d(4);
    Band5(index) = d(5);  
end

% X-M section
for KappaY=0+pi/L/RES:pi/L/RES:pi/L
    KappaY
    index = index + 1;
    KappaX = pi/L;
    
    W   = PWETE(KappaX,KappaY,N,L,r,epa,epb);
    
    d = eig(W);
    d = sort(sqrt(d));
    Band1(index) = d(1);
    Band2(index) = d(2);
    Band3(index) = d(3);
    Band4(index) = d(4);
    Band5(index) = d(5);
    
    WTM = PWETM(KappaX,KappaY,N,L,r,epa,epb);
    
    dTM = eig(WTM);
    dTM = sort(sqrt(dTM));
    Band1TM(index) = dTM(1);
    Band2TM(index) = dTM(2);
    Band3TM(index) = dTM(3);
    Band4TM(index) = dTM(4);
    Band5TM(index) = dTM(5);
    
end

% M-GAMMA section
for KappaX=pi/L-pi/L/RES:-pi/L/RES:0
    KappaX
    index = index + 1;
    KappaY = KappaX;
    
    W   = PWETE(KappaX,KappaY,N,L,r,epa,epb);
    
    d = eig(W);
    d = sort(sqrt(d));
    Band1(index) = d(1);
    Band2(index) = d(2);
    Band3(index) = d(3);
    Band4(index) = d(4);
    Band5(index) = d(5);
    
    WTM = PWETM(KappaX,KappaY,N,L,r,epa,epb);
    
    dTM = eig(WTM);
    dTM = sort(sqrt(dTM));
    Band1TM(index) = dTM(1);
    Band2TM(index) = dTM(2);
    Band3TM(index) = dTM(3);
    Band4TM(index) = dTM(4);
    Band5TM(index) = dTM(5);
    
end
    
MAX = max(max(Band3TM),max(Band3));
KAPPA_SCAN = zeros(size(Band1));
for i=1:2*RES+1
    KAPPA_SCAN(i) = (i-1)*pi/L/RES;
end

for j=i+1:3*RES+1
    KAPPA_SCAN(j) = (i-1)*pi/L/RES + (j-i)*pi/L/RES*sqrt(2);
end

KAPPA_SCAN;
plot(KAPPA_SCAN*L,Band1*L/2/pi);
hold;
plot(KAPPA_SCAN*L,Band2*L/2/pi);
plot(KAPPA_SCAN*L,Band3*L/2/pi);
plot(KAPPA_SCAN*L,Band4*L/2/pi);
plot(KAPPA_SCAN*L,Band5*L/2/pi);

plot(KAPPA_SCAN*L,Band1TM*L/2/pi,'--r');
plot(KAPPA_SCAN*L,Band2TM*L/2/pi,'--r');
plot(KAPPA_SCAN*L,Band3TM*L/2/pi,'--r');
plot(KAPPA_SCAN*L,Band4TM*L/2/pi,'--r');
plot(KAPPA_SCAN*L,Band5TM*L/2/pi,'--r');

plot(KAPPA_SCAN(RES+1)*L,0:MAX*L/2/pi/100:MAX*L/2/pi);
plot(KAPPA_SCAN(i)*L,0:MAX*L/2/pi/100:MAX*L/2/pi);
plot(KAPPA_SCAN(j)*L,0:MAX*L/2/pi/100:MAX*L/2/pi);

title('Blue-Solid Line: E-polarization, Red-Dashed Line: H-polarization');
text(1.5,0.05,['r/L=',num2str(r),', {\epsilon}_a=',num2str(epa),', {\epsilon}_b=',num2str(epb)]);


XLabel('{\kappa}_n','Fontsize',12);
YLabel('{\omega}_n','Fontsize',12);
hold off;

axis tight;



function W=PWETE(KappaX,KappaY,N,L,r,epa,epb);

for q=-N:N
    for p=-N:N
        for n=-N:N
            for m=-N:N
                S_MATRIX((p+N+1)+(q+N)*(2*N+1),(m+N+1)+(n+N)*(2*N+1)) = ...
                    SmnpqTE(m,n,p,q,KappaX,KappaY,L,r,epa,epb);
            end
        end
    end
end



W = S_MATRIX;



function Smnpq=SmnpqTE(m,n,p,q,KappaX,KappaY,L,r,epa,epb)

Smnpq = (((2*pi/L)^2)*(m^2+n^2)+(4*pi/L)*(m*KappaX+n*KappaY)+(KappaX^2+KappaY^2)) ...
    *Relative_Epsilon_Expansion(1/epa,1/epb,L,r,p-m,q-n);


function W=PWETM(KappaX,KappaY,N,L,r,epa,epb);

for q=-N:N
    for p=-N:N
        for n=-N:N
            for m=-N:N
                S_MATRIX((p+N+1)+(q+N)*(2*N+1),(m+N+1)+(n+N)*(2*N+1)) = ...
                    SmnpqTM(m,n,p,q,KappaX,KappaY,L,r,epa,epb);
            end
        end
    end
end

W = S_MATRIX;


function Smnpq=SmnpqTM(m,n,p,q,KappaX,KappaY,L,r,epa,epb)

Smnpq = (KappaX^2+KappaY^2+(2*pi/L)*((m+p)*KappaX+(n+q)*KappaY)+((2*pi/L)^2)*(m*p+q*n)) ...
    * Relative_Epsilon_Expansion(1/epa,1/epb,L,r,p-m,q-n);



function epmn=Relative_Epsilon_Expansion(epa,epb,L,r,m,n)

Gmn = (2*pi/L)*sqrt(m^2+n^2);

if r>0
    if Gmn==0 
        epmn = epb+(epa-epb)*pi*(r^2)/L^2;
    else
        epmn = epb*mydelta(r*Gmn) + (-epb+epa)*2*pi*(r^2)/(L^2)*besselj(1,Gmn*r)/(Gmn*r);
    end
else
    if Gmn==0
        epmn = epb;
    else
        epmn = 0;
    end
end

function output=mydelta(x)

output = 0;
if x==0
    output = 1;
end