% FDTD 2D with PML
% based on Sullivan 2D + PML
% and Lazzi 2D code
% Carlos J. Cela, March 2012

clear all;
close all;

e0 = 8.854e-12;
u0=4*pi*1.e-7;

% Model size (total, including PML)
ie=30;
je=30; 	

% Number of PML layers
npml=5;

%Lower corner of modeling area
ia=npml+1;
ja=npml+1; 	

%Upper corner of modeling area
ib=ie-npml;
jb=je-npml; 

% Number of materials
nmed=0; 

% Set the resolution, time step, frequency, omega, lamda, 
% and number of time steps 

ddx = 1.0e-3;
dt = ddx/(6.e8);

freq=1.9e9;
omega=2*pi*freq;
lamda=3.e8/freq;

% number of time steps per cycle
np=(2.*3.e8)/(freq*ddx);

% number of time steps for the simulation
nsteps =round(6.*np); 

% Source placement
ic = ie/2.;
jc = je/2.;

% Initialize the arrays

ez=zeros(ie,je);
dz=zeros(ie,je);
hx=zeros(ie,je);
hy=zeros(ie,je);
iz=zeros(ie,je);
ga=zeros(ie,je)+1.0;
gb=zeros(ie,je);

ihx=zeros(ie,je);
ihy=zeros(ie,je);

% Initialize the PML

fi1(1:ie)=0.;
gi2(1:ie)=1.;
fi2(1:ie)=1.;
gi3(1:ie)=1.;
fi3(1:ie)=1.;

fj1(1:je)=0.;
gj2(1:je)=1.;
fj2(1:je)=1.;
gj3(1:je)=1.;
fj3(1:je)=1.;

for i=1:npml
	xxn=1.*(npml+1-i)/(npml*1.);
	xn=.33*(xxn^3);
	gi2(i)=1./(1.+xn);
	gi2(ie-i+1)=1./(1.+xn);
	gi3(i)=(1.-xn)/(1.+xn);
	gi3(ie-i+1)=(1.-xn)/(1.+xn);

	xxn=(npml+1-i-0.5)/(npml*1.);
	xn=.33*(xxn^3);
	fi1(i)=xn;
	fi1(ie-i)=xn;
	fi2(i)=1./(1.+xn);
	fi2(ie-i)=1./(1.+xn);
	fi3(i)=(1.-xn)/(1.+xn);
	fi3(ie-i)=(1.-xn)/(1.+xn); 
end;

for j=1:npml
	xxn=1.*(npml+1-j)/(npml*1.);
	xn=.33*(xxn^3);
	gj2(j)=1./(1.+xn);
	gj2(je-j+1)=1./(1.+xn);
	gj3(j)=(1.-xn)/(1.+xn);
	gj3(je-j+1)=(1.-xn)/(1.+xn);
   
  xxn=(npml+1-j-0.5)/(npml*1.);
	xn=.33*(xxn^3);
	fj1(j)=xn;
	fj1(je-j)=xn;
	fj2(j)=1./(1.+xn);
	fj2(je-j)=1./(1.+xn);
	fj3(j)=(1.-xn)/(1.+xn);
	fj3(je-j)=(1.-xn)/(1.+xn);
end

% Define air
eps=1.;
cond=0.;
for i=1:ie,
  for j=1:je,
	  ga(i,j)=1./(eps+(cond*dt/e0));
    gb(i,j)=cond*dt/e0;
  end
end

% Start the simulation
figure()
hold on;
T = 0
for n=1:nsteps
  T = T+1;

  % Calculate the D field
  for i=2:ie-1,
    for j=2:je-1,
		  curl_h=(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
		  dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*.5*(curl_h);
    end
  end


  pulse=exp(-0.5*((T-80.)/20.)^2);
  dz(ic,jc)=pulse;

  % Calculate the E field
  for i=1:ie-1,
    for j=1:je-1,
      ez(i,j)=ga(i,j)*(dz(i,j)-iz(i,j));
	    iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j);
    end
  end

  % Calculate the H field
  for i=1:ie-1,
    for j=1:je-1,
	    curl_e=(-ez(i,j+1)+ez(i,j));
	    ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;
	    hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*.5*(curl_e+ihx(i,j));
    end
  end

  for i=1:ie-1,
    for j=1:je-1,
	    curl_e=(ez(i+1,j)-ez(i,j));
	    ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;
	    hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*.5*(curl_e+ihy(i,j));
    end
  end

  imagesc(ez, [0 0.5]);

  pause(0.001);
end





