FDTD_2_1

%%%%%%%% 1D FDTD simulation %%%%%%%% %%%%%%%% More on One-dimensional simulation 2_1 %%%%%%%% %%%%%%%% Reformulation using the flux density %%%%%%%%

KE=200;

ex=zeros(1,KE+1);

hy=zeros(1,KE+1);

ga=zeros(1,KE+1);

gb=zeros(1,KE+1);

dx=zeros(1,KE+1);

ix=zeros(1,KE+1);

kc=KE/2;

t0=40;

spread=12;

T=0;

NSTEP=300;

ex_low_m1=0;

ex_low_m2=0;

ex_high_m1=0;

ex_high_m2=0;

ddx=0.01; %set the cell size to 1cm

dt=ddx/(2*3e8); %calculate the time step

freq_in=7*1e8; %signal frequency

epsilon=4; %relative dielectric constant

kstart=100; %position of the dielectric medium

epsz=8.85419e-12; %absolute dielectric constant

sigma=0.04; %conductivity

%initialize to free space

for n=1:KE+1

ga(n)=1;

end

for n=kstart:KE+1

ga(n)=1/(epsilon+(sigma*dt/epsz));

gb(n)=sigma*dt/epsz;

end

for n=1:NSTEP

T=T+1;

%calculate dx field

for k=2:KE+1

dx(k)=dx(k)+0.5*(hy(k-1)-hy(k));

end

%insert source

pulse=exp(-0.5*(t0-T)^2/spread^2); %

sin(2*pi*freq_in*dt*T);

dx(5)=dx(5)+pulse; %%%%%%%%%%ex(50)=ex(50)+pulse与ex(50)=pulse的区别？？

%calculate E from dx

for k=2:KE+1

ex(k)=ga(k)*(dx(k)-ix(k));

ix(k)=ix(k)+gb(k)*ex(k);

end

%absorbing boundary condition

ex(1)=ex_low_m2;

ex_low_m2=ex_low_m1;

ex_low_m1=ex(2);

ex(KE) = ex_high_m2;

ex_high_m2 = ex_high_m1;

ex_high_m1 = ex(KE-1);

%calculate H field

for k=1:KE

hy(k)=hy(k)+0.5*(ex(k)-ex(k+1));

end

m=moviein(300);

i=1:KE+1;

subplot(2,1,1);

plot(i,ex);

axis([0 200 -1.5 1.5])

title(['电场' 't= ' num2str(n)]) ;

subplot(2,1,2);

plot(i,hy);

axis([0 200 -1.5 1.5])

title(['磁场''t= ' num2str(n)]) ;

m(n)=getframe(gcf);

%pause(0.005)

end