In [6]:
using Plots;
gadfly();
In [7]:
include("../fdtd/update.jl");
include("../fdtd/sources.jl");
include("../fdtd/boundaries.jl");
using update;
using sources;
In [8]:
#Global parameters
size = 200;
endTime = 500;
num_snaps = 250;
snap_step = div(endTime, num_snaps);
LOSS = 0.02;
eps1 = 9;
LossLayerStart = 180;
#Grid
# Magnetic
hy = zeros(size-1);
mu = ones(size-1);
chyh = ones(size);
chye = ones(size);
# Electric
ez = zeros(size);
eps = ones(size);
cezh = ones(size);
ceze = ones(size);
for i in 100:size
eps[i] = eps1;
end
for i in LossLayerStart:size
ceze[i] = (1.0 - LOSS) / (1.0 + LOSS);
cezh[i] = 1.0 / (1.0 + LOSS);
end
for i in LossLayerStart:size-1
chyh[i] = (1.0 - LOSS) / (1.0 + LOSS);
chye[i] = 1.0 / (1.0 + LOSS);
end
# output params
ez_snapshot = Array{Any}(num_snaps);
hy_snapshot = Array{Any}(num_snaps);
In [9]:
# Time steps
for time in 1:endTime
# Incident
# ez_inc, hy_inc = sources.gaussian_source(50, time);
delay = 30.
width = 100.
ez_inc = exp(-(time + 0.5 - (-0.5) - delay) * (time + 0.5 - (-0.5) - delay) / width);
hy_inc = exp(-(time - delay) * (time - delay) / width) / globals.imp0;
#
# Magnetic
#
# Interior update
update.update_magnetic_field!(ez, hy, mu, chyh, chye);
# TFSF
hy[49] -= hy_inc
#
# Electric
#
# ABC
# boundaries.trivial_abc!(ez, size);
boundaries.trivial_abc!(ez, 1, false);
# Interior update
update.update_electric_field!(ez, hy, eps, cezh, ceze);
# PEC
# boundaries.pec_boundary!(ez, size);
# TFSF
ez[50] += ez_inc
# Snapshots for animation
if mod(time, snap_step) == 0
ez_snapshot[div(time,snap_step)] = (time, copy(ez))
hy_snapshot[div(time,snap_step)] = (time, copy(hy).*globals.imp0)
end
end
In [11]:
anim = Animation()
for i = 1:num_snaps
p = plot(1:size, ez_snapshot[i][2], lab="Ez")
plot!(1:size-1, hy_snapshot[i][2], lab="Hy*imp0")
time = ez_snapshot[i][1]
plot!(ann=[(150, 1.5, "time =$time")])
plot!(ann=[(0, 1.1, "ABC")])
plot!(ann=[(80, 1.2, "Eps = 1")])
plot!(ann=[(100, 1.1, "Eps = 9")])
plot!([100, 100], [-2, 2])
plot!(ann=[(180, 1.1, "PML")])
plot!([180, 180], [-2, 2])
plot!(xlims=(1, 200), ylims=(-2, 2))
frame(anim, p)
end
gif(anim, "./07_impedance.gif", fps=15)
Out[11]:
In [ ]: