%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

lamb=1.55e-6;
Dlamb=1;%e-6;  %The lattice constant.
R0=1.12*Dlamb;  %中心缺陷

%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=61; %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=5; %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=2500; %Total number of Time Steps

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

R=0.48; %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.
  f=c/lamb;             
a=1;%e-6;   %The lattice constant.
W=W*(2*pi*c/a); %frequency
               
Dx=(2*Dlamb)/(NMlat-1); %Delta x.
%Ep_cell=ones(NMlat,NMlat)*e0;錯誤
Ep_cell=ones(NMlat,NMlat)*e0*ea;%  21*21
M=(NMlat-1)/2*Dx;
x=-M:Dx:M;%   21
N=sqrt(3/4)*M;
Dy=(N*2)/(NMlat-1); %Delta y.
Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c; %Time interval
y=-N:Dy:N;


[X,Y]=meshgrid(x,y);
X=X';%   21*21
Y=Y';%   21*21
flag=find(sqrt(X.^2+Y.^2)<R);%find elements' suffix of radius less than R.
Ep_cell(flag)=e0;%  21*21
flag1=find(sqrt((X+M).^2+Y.^2)<R);%find elements' suffix of radius less than R.left
Ep_cell(flag1)=e0;
flag2=find(sqrt((X-M).^2+Y.^2)<R);%find elements' suffix of radius less than R.right
Ep_cell(flag2)=e0;
flag3=find(sqrt((X+(M/2)).^2+(Y-N).^2)<R);%find elements' suffix of radius less than R.up1
Ep_cell(flag3)=e0;
flag4=find(sqrt((X-(M/2)).^2+(Y-N).^2)<R);%find elements' suffix of radius less than R.up2
Ep_cell(flag4)=e0;
flag5=find(sqrt((X+(M/2)).^2+(Y+N).^2)<R);%find elements' suffix of radius less than R.down1
Ep_cell(flag5)=e0;
flag6=find(sqrt((X-(M/2)).^2+(Y+N).^2)<R);%find elements' suffix of radius less than R.down2
Ep_cell(flag6)=e0;
%surf(X,Y,Ep_cell/e0);
%shading interp;
% view(0,90);
% axis('square');
% axis off;
 Ep=repmat(Ep_cell,MLatx,MLaty);% copy Ep_cell in 11*11 lattice matrix.


%toc
 if WantToSeeEp==1
   x=0:Dx:(NMlat*MLatx-1)*Dx;% 231 Dx start from 0.
   x=x-(NMlat*MLatx-1)*Dx/2;% value from -(NMlat*MLatx-1)*Dx/2 to (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';
   flag=find(sqrt(X.^2+Y.^2)<R0);%find elements' suffix of radius less than R.
   Ep(flag)=e0;%  21*21
   surf(X,Y,Ep/e0);
   shading interp;
   colorbar     %display a colorful figure.
   %colormap([1 1 1;0 0 0]) %disply a black and white figure.
   view(0,90);
   axis([min(x), max(x),min(y), max(y)])
   axis('square');  %adjust a figure's display 
   axis off;
   disp('Press any key to continue...');
   pause
end%End of defining the Ep.