%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2D TM_FDTD
%                    by M.Takahashi     2003.6.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

Nx=30, Ny=30
NStep=5, NStop=100

delx=1./20., dely=1./20.    %unit='m'
eps0=8.854E-12, xmu0=1.2566306E-6

C=1.0/sqrt(xmu0*eps0)

dTXI=C/delx, dTYI=C/dely
dT=1./sqrt(dTXI^2+dTYI^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Incidet Waveform Parameters
Beta=32.0
BetadT =  Beta*dT, Alpha=(1./(BetadT/4.0))^2, Priod = 2.0*BetadT
Ife=15; Jfe=15    %feeding position
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Graphic Position
scrsz = get(0,'ScreenSize');
figure(1), clf, set(1,'Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
figure(2), clf, set(2, 'Position',[scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
M = moviein(fix(NStop/NStep)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Initialize
Ez=zeros(Nx,Ny), Hx=zeros(Nx-1,Ny-1), Hy=zeros(Nx-1,Ny-1), Ez1=zeros(Nx,Ny), Ez2=zeros(Nx,Ny)
CEz=1., CEzlx=dT/(eps0*delx), CEzly=dT/(eps0*dely), CHylx=dT/(xmu0*delx), CHxly=dT/(xmu0*dely)
%   Coefficient for Mur's 2nd ABC
Cxd=(C*dT-delx)/(C*dT+delx), Cyd=(C*dT-dely)/(C*dT+dely)
Cxu=Cxd, Cyu=Cyd, Cxx=2.*delx/(C*dT+delx), Cyy=2.*dely/(C*dT+dely)
Cxfyd=delx*C*dT*C*dT/(2.*dely*dely*(C*dT+delx)), Cyfxd=dely*C*dT*C*dT/(2.*delx*delx*(C*dT+dely))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%-------MAIN LOOP---------------------------------------------
T=0.0
Nf=0

for N=1: NStop
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Incident Source
	Source(N)=0.
	if (T >= 0.) & (T<= Priod)
		Source(N)=exp(-Alpha*((T-BetadT)^2))
 	end    
 	Tau(N)=T
    
    figure(1), plot(Tau,Source, 'o'), axis([0 dT*NStop 0 1]), xlabel('Time (s)'), ylabel('Source Amplitude')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ez Field
%  	for J=2: Ny-1
%    	for I=2: Nx-1
%        	Ez(I,J)= CEz*Ez(I,J) + CEzlx*(Hy(I,J)-Hy(I-1,J)) - CEzly*(Hx(I,J)-Hx(I,J-1))
%    	end
%  	end
   	Ez(2:Nx-1,2:Ny-1)= CEz*Ez(2:Nx-1,2:Ny-1) + CEzlx*(Hy(2:Nx-1,2:Ny-1)-Hy(1:Nx-2,2:Ny-1)) - CEzly*(Hx(2:Nx-1,2:Ny-1)-Hx(2:Nx-1,1:Ny-2))
   
   	Ez(Ife,Jfe)=Source(N)    % Source excitation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot
	if mod(N,NStep)==0
      	Nf=Nf+1
      	figure(2), surf(Ez), caxis([0 0.2]), colorbar, axis([0 Nx 0 Ny -0.3 0.5])
        xlabel('Nx'), ylabel('Ny'), zlabel('Ez')
%		shading interp
		M(:,Nf) = getframe; 
	end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mur's 2nd ABC
    Ez(1,2)=Ez1(2,2)+Cxd*(Ez(2,2)-Ez1(1,2))
    Ez(Nx,2)=Ez1(Nx-1,2)+Cxu*(Ez(Nx-1,2)-Ez1(Nx,2))
    Ez(1,Ny-1)=Ez1(2,Ny-1)+Cxd*(Ez(2,Ny-1)-Ez1(1,Ny-1))
    Ez(Nx,Ny-1)=Ez1(Nx-1,Ny-1)+Cxu*(Ez(Nx-1,Ny-1)-Ez1(Nx,Ny-1))

    Ez(1,3:Ny-2)=-Ez2(2,3:Ny-2)+Cxd*(Ez(2,3:Ny-2)+Ez2(1,3:Ny-2)) +Cxx*(Ez1(1,3:Ny-2)+Ez1(2,3:Ny-2)) +Cxfyd*(Ez1(1,4:Ny-1)-2.*Ez1(1,3:Ny-2)+Ez1(1,2:Ny-3) +Ez1(2,4:Ny-1)-2.*Ez1(2,3:Ny-2)+Ez1(2,2:Ny-3))
    Ez(Nx,3:Ny-2)=-Ez2(Nx-1,3:Ny-2)+Cxd*(Ez(Nx-1,3:Ny-2)+Ez2(Nx,3:Ny-2)) +Cxx*(Ez1(Nx,3:Ny-2)+Ez1(Nx-1,3:Ny-2)) +Cxfyd*(Ez1(Nx,4:Ny-1)-2.*Ez1(Nx,3:Ny-2)+Ez1(Nx,2:Ny-3) +Ez1(Nx-1,4:Ny-1)-2.*Ez1(Nx-1,3:Ny-2)+Ez1(Nx-1,2:Ny-3))
    
    Ez(2,1)=Ez1(2,2)+Cyd*(Ez(2,2)-Ez1(2,1))
    Ez(2,Ny)=Ez1(2,Ny-1)+Cyu*(Ez(2,Ny-1)-Ez1(2,Ny))
    Ez(Nx-1,1)=Ez1(Nx-1,2)+Cyd*(Ez(Nx-1,2)-Ez1(Nx-1,1))
    Ez(Nx-1,Ny)=Ez1(Nx-1,Ny-1)+Cyu*(Ez(Nx-1,Ny-1)-Ez1(Nx-1,Ny))
    
    Ez(3:Nx-2,1)=-Ez2(3:Nx-2,2)+Cyd*(Ez(3:Nx-2,2)+Ez2(3:Nx-2,1)) +Cyy*(Ez1(3:Nx-2,1)+Ez1(3:Nx-2,2)) +Cyfxd*(Ez1(4:Nx-1,1)-2.*Ez1(3:Nx-2,1)+Ez1(2:Nx-3,1) +Ez1(4:Nx-1,2)-2.*Ez1(3:Nx-2,2)+Ez1(2:Nx-3,2))
    Ez(3:Nx-2,Ny)=-Ez2(3:Nx-2,Ny-1)+Cyd*(Ez(3:Nx-2,Ny-1)+Ez2(3:Nx-2,Ny)) +Cyy*(Ez1(3:Nx-2,Ny)+Ez1(3:Nx-2,Ny-1)) +Cyfxd*(Ez1(4:Nx-1,Ny)-2.*Ez1(3:Nx-2,Ny)+Ez1(2:Nx-3,Ny) +Ez1(4:Nx-1,Ny-1)-2.*Ez1(3:Nx-2,Ny-1)+Ez1(2:Nx-3,Ny-1))

    Ez2= Ez1    % Stored Ez Value at T=(n-2)*dt
    Ez1= Ez     % Stored Ez Value at T=(n-1)*dt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	T=T+dT/2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hx Field
%  	for J=1: Ny-1
%    	for I=2: Nx-1
%      		Hx(I,J)=Hx(I,J) - CHxly*(Ez(I,J+1)-Ez(I,J))
%    	end
%  	end
	Hx(2:Nx-1,1:Ny-1)=Hx(2:Nx-1,1:Ny-1) - CHxly*(Ez(2:Nx-1,2:Ny)-Ez(2:Nx-1,1:Ny-1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hy Field
%  	for J=2: Ny-1
%    	for I=1: Nx-1
%      		Hy(I,J)=Hy(I,J) + CHylx*(Ez(I+1,J)-Ez(I,J))
%    	end
%  	end
    Hy(1:Nx-1,2:Ny-1)=Hy(1:Nx-1,2:Ny-1) + CHylx*(Ez(2:Nx,2:Ny-1)-Ez(1:Nx-1,2:Ny-1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	T=T+dT/2

end

movie(M,3)
