fd1d2.m

2D FDTD Simulation program implemented in MATLAB

% Fd1d2.m Dielectric slab

clear all

NN = 200;
kc = NN/2;
k0 = 100;

epsz = 8.85419*10^-12;
ddx = .01;
dt = ddx/(2*3*10^8);

t0 = 50;
spread = 15;

ex = zeros(1,NN);
hy = zeros(1,NN);
ca = zeros(1,NN);
cb = zeros(1,NN);
cond = zeros(1,NN);
barr = zeros(1,NN);

elm1 = 0.;
elm2 = 0.;
ehm1 = 0.;
ehm2 = 0.;

% Initialize
for k=1:NN
ca(k) = 1.;
cb(k) = 0.5;
end

epsilon = 5;
sigma = 0.01;
eaf  = dt*sigma/(2*epsz*epsilon);

for k=80:120
cond(k) = sigma;
ca(k) = (1 -eaf)/(1+eaf);
cb(k) = .5/(epsilon*(1+eaf));
barr(k) = .5;
end

freq = zeros(2,1);
real_pt = zeros(2,NN);
imag_pt = zeros(2,NN);
amp = zeros(2,NN);

freq(1) = 50*1e6;
freq(2) = 500*1e6;

T = 0;
n_step = 1;
while n_step > 0
n_step = input('How many time steps -->');

for n=1:n_step
T = T + 1;

ex(1) = elm2;
elm2 = elm1
elm1 = ex(2)
for k=2:NN
ex(k) = ca(k)*ex(k) + cb(k)*( hy(k-1) - hy(k) );
end
ex(NN) = ehm2;
ehm2 = ehm1;
ehm1 = ex(NN-1)

for k=2:NN
for n=1:2
real_pt(n,k) = real_pt(n,k) + cos(2*pi*T*dt*freq(n))*ex(k);
imag_pt(n,k) = imag_pt(n,k) + sin(2*pi*T*dt*freq(n))*ex(k);
end
end

pulse = exp(-.5*((t0-T)/spread)^2 );
%pulse = sin(2*pi*T/50);
pulse
if T < 120
ex(2) = pulse;
end

for k=1:NN-1
hy(k) =  hy(k) + .5*( ex(k) - ex(k+1) );
end

end

subplot(3,1,1)
plot(ex)
hold on
plot(barr,'m--')
hold off
axis( [ 1 NN -.5 1.1 ])
TT = text(150, .7 , ['T = ',num2str(T)]);
set(TT,'fontsize',12)
TT = text(87, .3 , [' eps = ',num2str(epsilon)]);
set(TT,'fontsize',12)
TT = text(87, -.3 , ['sig = ',num2str(sigma)]);
set(TT,'fontsize',12)
set(gca,'fontsize',12)
ylabel('E field')

for k=1:NN
for n=1:2
amp(n,k) = cond(k)*(real_pt(n,k)^2 + imag_pt(n,k)^2);
end
end

subplot(6,1,3)
plot(amp(1,:))
set(gca,'fontsize',12)
TT = text(130, 3, ['50 MHz']);
set(TT,'fontsize',12)
axis( [ 1 NN 0 7 ])
set(gca,'XTick',[])
ylabel('SAR')

subplot(6,1,4)
plot(amp(2,:))
set(gca,'fontsize',12)
TT = text(130, 3, ['500 MHz']);
axis( [ 1 NN 0 7 ])
set(TT,'fontsize',12)
ylabel('SAR')

saveas(gcf,'pulse.bmp')

end