%This program calculate the waveguide by a photonic crystal.
%For two dimension case.
%Only for TM case(Ez,Hx,Hy), in this version試驗版本.
clear;
%tic
%Initial parameters and other things.
W=0.36; %Normalized frequency 就是能帶圖中的縱坐標，單位為（a/lamda）

%The following parameters are control parameters.
WaveGuide=1; %If this program is for Wave Guide? If so, please specify 1.
IsMovie=0; %If you want to play movie, please use 1.
IsFigure=1; %If it will plot the figures? If so, please specify 1.
WantToSeeEp=1; %Do you want to see the distrubution of Ep（TE極化，S極化）? If so, please specify 1.
%End of defining control parameters

MLatx=11; %How many Lattice cell in x direction.光子晶體柵格
MLaty=11; %How many Lattice cell in y direction.

NMlat=21; %The gird number in each Lattice Cell. 每個柵格中包含的fdtd網格個數
          %SHOULD BE ODD INTEGER!!!!奇數整數個,目的就是為了將光子晶體的原點定位在網格中心
if mod(NMlat,2)==0
   NMlat=NMlat+1;
end       %Force it to be a odd integer!

NTx=MLatx*NMlat+1; %It is the number of the Grid along x axis.總的fdtd網格坐標點（所以要加1，呵呵）
NTy=MLaty*NMlat+1; %It is the number of the Grid along y axis.

if WaveGuide==1
   Nrow=4; %The row number of coloumns between the PML boundary and the waveguide.(此處定義，現在還不知所措)
end

NPML=12; %How many PML layers will be used in our computation.PML的層數

NTimeSteps=2500; %Total number of Time Steps

Meach=20; %Define the interval for plot figures if IsFigure==1.
           %This also works for saving intervals.每隔20個時間步顯示一次。
           

R=0.2; %The radius of dielectric columns，光子晶體的半徑，可調
ea=11.4; %The dielectric constant of these columns.

Zmax=0.6;    %The maximum value for z axis when plotting figures.
Colormax=0.6; %The maximum value for colormap when plotting figures.


%Some constants
mu0=4*pi*1.0e-7; %Epsilon Zero, if using Gauss Unit, it equals to 1.
e0=8.85*1e-12;   %Mu Zero, if using Gauss Unit, it equals to 1.
c=1/sqrt(mu0*e0); %The light speed.
factor=mu0/e0; %The factor between conductivity and permeability. 
               %Permeability=Conductivity*factor, in PML.PML中的匹配條件
               
a=1;%e-6;   %The lattice constant.柵格常數。
W=W*(2*pi*c/a); %frequency
               
Dx=a/NMlat; %Delta x.dx＝1/(一個柵格中的網格數)，此程序是先給光子晶體的數目，再給fdtd網格分辨率。
Dy=Dx; %Delta y.
Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c; %Time interval

%tic
%In the following partm we define the dielectric
%constants:，定義FDTD中光子晶體與空氣的介電常數
Ep=ones(NTx-1,NTy-1)*e0;  %FDTD總的網格中的介質常數，真空。
Ep_cell=ones(NMlat,NMlat)*e0;%每個柵格中的fdtd介電常數，真空。

x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;%先對柵格操作，然后進行repmat，復制并平鋪陣列。
[X,Y]=meshgrid(x);%把x方向上的向量拓展為X,Y兩個矩陣。
X=X';
Y=Y';
flag=find(sqrt(X.^2+Y.^2)<R); %尋找光子晶體，并賦給介電常數
Ep_cell(flag)=e0*ea;  
Ep=repmat(Ep_cell,MLatx,MLaty);%拓展到真個空間了。

if WaveGuide==1
   eb=1; %The dielectric constant in the waveguide
   Ep(1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*eb;
   Ep((MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:MLaty*NMlat)=e0*eb;
   Ep((MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=Ep_cell;
   Ep((MLatx-Nrow-2)*NMlat+1:(MLatx-Nrow-1)*NMlat,(Nrow+1)*NMlat+1:(Nrow+2)*NMlat)=e0*eb;
end
%toc

if WantToSeeEp==1
   x=0:Dx:(NMlat*MLatx-1)*Dx;
   x=x-(NMlat*MLatx-1)*Dx/2;
   y=0:Dy:(NMlat*MLaty-1)*Dy;
   y=y-(NMlat*MLaty-1)*Dy/2;
   [X,Y]=meshgrid(x,y);
   X=X';
   Y=Y';
   surf(X,Y,Ep/e0);
   shading interp;
   view(0,90);
   axis([min(x), max(x),min(y), max(y)])
   axis off;
   disp('Press any key to continue...');
   pause
end
%End of defining the Ep.


if IsFigure==1
   %Define the X Y coordinate for figures
   x=0:Dx:(NMlat*MLatx-1+2*NPML)*Dx;
   x=x-(NMlat*MLatx-1+2*NPML)*Dx/2;
   y=0:Dy:(NMlat*MLaty-1+2*NPML)*Dy;
   y=y-(NMlat*MLaty-1+2*NPML)*Dy/2;
   [X,Y]=meshgrid(x,y);
   X=X';
   Y=Y';
end

%Define the Ez, Hx, Hy, which are in the inside region.
Ez=zeros(NTx-1,NTy-1); 
Hx=zeros(NTx-1,NTy);
Hy=zeros(NTx,NTy-1);

%Parameters about PML:
n=4; %The order of the polynomial that decribes the conductivity profile.
R=1e-10;
Delta=NPML*Dx;
SigmaMax=-(n+1)*e0*c*log(R)/(Delta*2); 
NUM=NPML*2:-1:1;

Sigmax=SigmaMax*((NUM*Dx/2+Dx/2).^(n+1)-(NUM*Dx/2-Dx/2).^(n+1))/(Delta^n*Dx*(n+1));
Sigmay=Sigmax;
SigmaBound=SigmaMax*(Dx/2).^(n+1)/(Delta^n*Dx*(n+1));
%Sigmax=SigmaMax*(NUM/6).^(n+1); %Another way, the definition.

EzxPML1=zeros(NPML,NPML);
EzyPML1=zeros(NPML,NPML);
HxPML1=zeros(NPML,NPML);
HyPML1=zeros(NPML,NPML); %Zone 1

Sigmax_z1=repmat(Sigmax(2:2:NPML*2)',1,NPML);
Sigmax_x1=repmat(Sigmax(2:2:NPML*2)',1,NPML);
Sigmax_y1=repmat(Sigmax(1:2:NPML*2-1)',1,NPML); 
Sigmay_z1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1));
Sigmay_x1=fliplr(repmat(Sigmax(1:2:NPML*2-1),NPML,1));
Sigmay_y1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1)); %Zone 1

EzxPML2=zeros(NPML,NPML);
EzyPML2=zeros(NPML,NPML);
HxPML2=zeros(NPML,NPML);
HyPML2=zeros(NPML,NPML); %Zone 2

Sigmax_z2=flipud(Sigmax_z1);
Sigmax_x2=flipud(Sigmax_x1);
Sigmax_y2=flipud(Sigmax_y1); 
Sigmay_z2=Sigmay_z1;
Sigmay_x2=Sigmay_x1;
Sigmay_y2=Sigmay_y1; %Zone 2

EzxPML3=zeros(NPML,NPML);
EzyPML3=zeros(NPML,NPML);
HxPML3=zeros(NPML,NPML);
HyPML3=zeros(NPML,NPML); %Zone 3

Sigmax_z3=Sigmax_z1;
Sigmax_x3=Sigmax_x1;
Sigmax_y3=Sigmax_y1;
Sigmay_z3=fliplr(Sigmay_z1);
Sigmay_x3=fliplr(Sigmay_x1);
Sigmay_y3=fliplr(Sigmay_y1);  %Zone 3

EzxPML4=zeros(NPML,NPML);
EzyPML4=zeros(NPML,NPML);
HxPML4=zeros(NPML,NPML);
HyPML4=zeros(NPML,NPML); %Zone 4

Sigmax_z4=flipud(Sigmax_z1);
Sigmax_x4=flipud(Sigmax_x1);
Sigmax_y4=flipud(Sigmax_y1); 
Sigmay_z4=fliplr(Sigmay_z1);
Sigmay_x4=fliplr(Sigmay_x1);
Sigmay_y4=fliplr(Sigmay_y1); %Zone 4

EzxPMLA=zeros(NTx-1,NPML);
EzyPMLA=zeros(NTx-1,NPML);
HxPMLA=zeros(NTx-1,NPML);
HyPMLA=zeros(NTx,NPML); %Zone A

Sigmay_zA=repmat(Sigmay_z1(1,:),NTx-1,1);
Sigmay_xA=repmat(Sigmay_x1(1,:),NTx-1,1);
Sigmay_yA=repmat(Sigmay_y1(1,:),NTx,1); %Zone A

EzxPMLB=zeros(NTx-1,NPML);
EzyPMLB=zeros(NTx-1,NPML);
HxPMLB=zeros(NTx-1,NPML);
HyPMLB=zeros(NTx,NPML); %Zone B

Sigmay_zB=fliplr(Sigmay_zA);
Sigmay_xB=fliplr(Sigmay_xA);
Sigmay_yB=fliplr(Sigmay_yA); %Zone B

EzxPMLC=zeros(NPML,NTy-1);
EzyPMLC=zeros(NPML,NTy-1);
HxPMLC=zeros(NPML,NTy);
HyPMLC=zeros(NPML,NTy-1);%Zone C

Sigmax_zC=repmat(flipud(Sigmay_z1(1,:)'),1,NTy-1);
Sigmax_xC=repmat(flipud(Sigmay_x1(1,:)'),1,NTy);
Sigmax_yC=repmat(flipud(Sigmay_y1(1,:)'),1,NTy-1); %Zone C

EzxPMLD=zeros(NPML,NTy-1);
EzyPMLD=zeros(NPML,NTy-1);
HxPMLD=zeros(NPML,NTy);
HyPMLD=zeros(NPML,NTy-1);%Zone D

Sigmax_zD=flipud(Sigmax_zC);
Sigmax_xD=flipud(Sigmax_xC);
Sigmax_yD=flipud(Sigmax_yC); %Zone D

if IsMovie==1
   Movie=moviein(NTimeSteps/Meach+1);
   Mnum=1;
   EzAll=[EzxPML3+EzyPML3,EzxPMLC+EzyPMLC,EzxPML1+EzyPML1;
          EzxPMLB+EzyPMLB,Ez,EzxPMLA+EzyPMLA;
          EzxPML4+EzyPML4,EzxPMLD+EzyPMLD,EzxPML2+EzyPML2];
   HxAll=[HxPML3,HxPMLC,HxPML1;HxPMLB,Hx,HxPMLA;HxPML4,HxPMLD,HxPML2];
   HyAll=[HyPML3,HyPMLC,HyPML1;HyPMLB,Hy,HyPMLA;HyPML4,HyPMLD,HyPML2];
   
   figure(1)
   surf(X,Y,real(EzAll))
   shading interp;
   axis([-0.5*a*MLatx-Dx*NPML 0.5*a*MLatx+Dx*NPML -0.5*a*MLaty-Dx*NPML 0.5*a*MLaty+Dx*NPML -Zmax Zmax])
   zlabel('Ez');
   Movie(:,1)=getframe;
end

if WaveGuide==1
   Npy=NMlat*Nrow+round((NMlat+1)/2);
   Npx=NMlat*Nrow+round((NMlat+1)/2);
else
   Npx=round(NTx/2);
   Npy=round(NTy/2);
end

expboundary=exp(-SigmaBound*Dt/e0);

tic
for m=1:NTimeSteps
   %tic
   %The following part is for the inside region. 
   %H components.
   Hx(:,2:NTy-1)=Hx(:,2:NTy-1)-Dt*(Ez(:,2:NTy-1)-Ez(:,1:NTy-2))/(Dy*mu0);
   Hy(2:NTx-1,:)=Hy(2:NTx-1,:)+Dt*(Ez(2:NTx-1,:)-Ez(1:NTx-2,:))/(Dx*mu0);
   
   Hx(:,NTy)=expboundary*Hx(:,NTy)-(1-expboundary)*...
      (EzxPMLA(:,1)+EzyPMLA(:,1)-Ez(:,NTy-1))/(SigmaBound*factor*Dy);     %Boundary A
   Hx(:,1)=expboundary*Hx(:,1)-(1-expboundary)*...
      (Ez(:,1)-EzxPMLB(:,NPML)-EzyPMLB(:,NPML))/(SigmaBound*factor*Dy);  %Boundary B

   Hy(1,:)=expboundary*Hy(1,:)+(1-expboundary)*...