sqideal.m

用平面波展开法模拟2维光子晶体的能带结构,FDTD法模拟2维光子晶体能量传输

clear all;
warning off;
epsa=8.9;
epsb=1;
a=1.0;
R=0.2*a;
i=sqrt(-1);
f=pi*R^2/a^2;
NumberofCell=1;
a1=a;
a2=a*i;
b1=2*pi/a/NumberofCell;
b2=2*pi/a/NumberofCell*i;
n=4;
display('fourier transform-----');
tic;
NumberofPW=(2*n+1)^2;
mind=(-2*n:2*n)'+2*n+1;
%aaaa=ones(1,size(mind))
%aaa=mind(:,ones(1,size(mind)));
mind=mind(:,ones(1,size(mind)))-2*n-1;
GG=mind'*b1+mind*b2;
eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R);
eps21(2*n+1,2*n+1)=epsb+f*(epsa-epsb);
zz=[-2 -2;-2 -1;-2 0;-2 1;-2 2;...
    -1 -2;-1 -1;-1 0;-1 1;-1 2;
    0 -2;0 -1;       0 1;0,2;
    1 -2;1 -1;1 0;1 1;1 2;
    2,-2;2,-1;2 0;2,1;2,2]*[a1 a2].';
eps20=zeros(length(eps21));
for x=1:length(zz)
    eps20=eps20+exp(i*(real(GG).*real(zz(x))+imag(GG).*imag(zz(x)) )).*eps21;
end
ff=pi*R.^2*length(zz)/(NumberofCell^2*a^2);
eps20=eps20./NumberofCell^2;
eps20(2*n+1,2*n+1)=epsb+ff*(epsa-epsb);
count=1;
for y=-n:n;
    for x=-n:n;
        G(count)=x*b1+y*b2;
        r(count,:)=[x,y];
        count=count+1;
    end
end

display('building ');
for x=1:NumberofPW;
    for y=x:NumberofPW;
        b=r(x,:)-r(y,:)+2*n+1;
        eps2(x,y)=eps20(b(2),b(1));
        eps2(y,x)=eps2(x,y);
    end
end

k1=2*pi/a*0.5.*(0:0.1:1);
k2=2*pi/a*(0.5+(0.1:0.1:1).*0.5*i);
k3=2*pi/a*(0.5+0.5*i).*(0.9:-0.1:0);
k0=[k1 k2 k3];

display('now calculate the eigen value');
eps2=inv(eps(2));
if max(max(real(eps2)))>10^7*max(max(imag(eps2)))
    display('symmetric');
    esp2=real(eps2);
else
    display('imag part is not zero');
    stop;
end 
option.disp=0;
counter=1;
for ii=1:length(k0);
    k=k0(ii);
    M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*eps2;
    E=sort(abs(eig(M)));
    freq(:,counter)=sqrt(abs(E(1:10))).*a./2./pi;
    display(sprintf('calculation of k=%f+%f*i is defined',real(k),imag(k)));
    counter=counter+1;
end
toc;
[max(freq(1,:)),min(freq(2,:))]

tmpx=1:length(k0);
plot(tmpx,freq,'o-','linewidth',2)
grid on;
axis([1 31 0 1]);