fdtd 1d

%
% ADI FDTD 1D test
% one dimension wave wave propagate toward x and polarize in z direction(
% Ez Hy )
% by sai1ormoon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the adi algorithm is 
%   
%  t->t+1/2
%  update Ez (implicit) 
%  update Hy (explicit) 
%  set periodic BC
%  
%  t+1/2->t+1
%  update Ez  (explicit) 
%  add source
%  update Hy (explicit) 
%  set periodic BC
% 
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
%grids
Nx=100;
%Ez field and Hy field
Ez=zeros(1,Nx+1);
Hy=zeros(1,Nx);


%phisical parameter
cc=2.99792458e8;
muz0=4*pi*1e-7;
epsz0=1/cc/cc/muz0;

%source parameter  lambda and frequency 
lambda = 500e-9;
freq = cc/lambda;
wavenum=freq*2*pi/cc;

%grid size
dx=lambda/20.0;
%courant condition
courant=0.8; 
dt = dx/cc/courant;


%%for Ez implicit in free space
Aez=zeros(1,Nx+1);
Bez=zeros(1,Nx+1);
Cez=zeros(1,Nx+1);
G=zeros(1,Nx+1);
P=zeros(1,Nx+1);

for(p=1:Nx+1)
    Aez(p)= -dt^2/4/epsz0/muz0/dx^2;
    Bez(p)= 1+2*dt^2/4/epsz0/muz0/dx^2;
    Cez(p)= -dt^2/4/epsz0/muz0/dx^2;
end

%fdtd main loop
for(t=1:600)

    %update Ez t->t+1/2 (implicit)
    for(e=2:Nx)
    P(e)=Ez(e)+dt/2/epsz0/dx*(Hy(e)-Hy(e-1));
    end
    %solve  A*Ez = P     A is tri-diagonal matrix  ,Ez is field , the  solver
    %algorithm is from 'Numerical Recipes to solve a tridiagonal matrix'
    bet = Bez(2);
    Ez(2)=P(2)/bet;
    for(i = 3:Nx)
        G(i)=Cez(i-1)/bet;
        bet=Bez(i)-Aez(i)*G(i); 
        Ez(i)=(P(i)-Aez(i)*Ez(i-1))/bet;
    end 
    for(i=Nx-1:-1:2)
        Ez(i)=Ez(i)-G(i+1)*Ez(i+1);
    end

    %update Hy t->t+1/2  (explicit)
    for(e=1:(Nx-1))
     Hy(e)=Hy(e)+dt/2/muz0/dx*(Ez(e+1)-Ez(e));
    end

    %simplify boundary using periodic boundary condition
    Ez(1)=Ez(Nx);
    Hy(Nx)=Hy(1);    
    
    %update Ez  t+1/2->t+1
    
    for(e=2:Nx)
     Ez(e)=Ez(e)+dt/2/epsz0/dx*(Hy(e)-Hy(e-1));
    end

    %source point 
    Ez(50)=Ez(50)+exp(-((t-50)/13)^2);
    %update Hy  t+1/2->t+1
    for(e=1:(Nx-1))
    Hy(e)=Hy(e)+dt/2/muz0/dx*(Ez(e+1)-Ez(e));
    end

    %periodic boundary
    Ez(1)=Ez(Nx);
    Hy(Nx)=Hy(1);

    %plot
    if mod(t,5)==0;

    rtime=num2str(round(t*dt/1.0e-15));

    subplot(2,1,1),plot(Ez),axis([0  (Nx) -1 1 ] );
    title(['time = ',rtime,' fs']);
    ylabel('EZ');

    subplot(2,1,2);
    plot(Hy);
    axis([ 0 Nx -3.0e-3 3.0e-3]);
    title(['time = ',rtime,' ns']);
    xlabel('x (meters)');
    ylabel('HY');
    pause(0.1);
    end

end