%This program calculate the wave guide by a photonic crystal.
%For two dimension case.
%Only for TM case, in this version.

%Trial Version

clear;
%tic
%Initial parameters and other things.

W=0.36; %Normalized frequency


%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? 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. 
          %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.
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.

NTimeSteps=1000; %Total number of Time Steps

Meach=20; %Define the interval for plot figures if IsFigure==1.
           %This also works for saving intervals.
           

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.
               
a=1;%e-6;   %The lattice constant.
W=W*(2*pi*c/a); %frequency
               
Dx=a/NMlat; %Delta x.
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:
Ep=ones(NTx-1,NTy-1)*e0; 

Ep_cell=ones(NMlat,NMlat)*e0;
x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;
[X,Y]=meshgrid(x);
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