In [1]:
using Plots;
gadfly();
In [2]:
include("../fdtd/update.jl");
include("../fdtd/sources.jl");
include("../fdtd/boundaries.jl");
using update;
using sources;
In [3]:
#Global parameters
size = 1000;
endTime = 10000;
num_snaps = 200;
snap_step = div(endTime, num_snaps);
# Incident
inc_pos = 200;
# Material
n1 = 1;
n3 = 4;
n2 = sqrt(n3*n1);
eps1 = n1^2;
eps2 = n2^2;
eps3 = n3^2;
wavelength = 800
q_wavelength = div(wavelength / n2, 4)
#Grid
# Magnetic
hy = zeros(size-1);
hy_0 = zeros(size-1);
mu = ones(size-1);
chyh = ones(size);
chye = ones(size);
# Electric
ez = zeros(size);
ez_0 = zeros(size);
eps = ones(size) * eps1;
eps_0 = ones(size) * eps1;
cezh = ones(size);
ceze = ones(size);
#for i in 110:170
# eps[i] = eps1;
#end
for i in div(size, 2):div(size, 2)+(q_wavelength)
eps[i] = eps2;
end
for i in div(size, 2)+(q_wavelength):size
eps[i] = eps3;
end
rightBound = boundaries.setup_first_order_abc(eps, mu, size, true)
leftBound = boundaries.setup_first_order_abc(eps, mu, 1, false)
rightBound_0 = boundaries.setup_first_order_abc(eps_0, mu, size, true)
leftBound_0 = boundaries.setup_first_order_abc(eps_0, mu, 1, false)
# output params
ez_snapshot = Array{Any}(num_snaps);
hy_snapshot = Array{Any}(num_snaps);
ez_0_snapshot = Array{Any}(num_snaps);
hy_0_snapshot = Array{Any}(num_snaps);
In [4]:
# Time steps
for time in 1:endTime
# Incident
ez_inc = sin(2*pi/(wavelength) * (time-1))*exp(-1000/time)
#
# Magnetic
#
# Interior update
update.update_magnetic_field!(ez, hy, mu, chyh, chye);
update.update_magnetic_field!(ez_0, hy_0, mu, chyh, chye);
#
# Electric
#
# Interior update
update.update_electric_field!(ez, hy, eps, cezh, ceze);
update.update_electric_field!(ez_0, hy_0, eps_0, cezh, ceze);
# ABC
boundaries.first_order_diff_abc!(ez, leftBound)
boundaries.first_order_diff_abc!(ez, rightBound)
boundaries.first_order_diff_abc!(ez_0, leftBound_0)
boundaries.first_order_diff_abc!(ez_0, rightBound_0)
# Incident
ez[inc_pos] += ez_inc / sqrt( eps[inc_pos] * mu[inc_pos])
ez_0[inc_pos] += ez_inc / sqrt( eps[inc_pos] * mu[inc_pos])
# 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)
ez_0_snapshot[div(time,snap_step)] = (time, copy(ez_0))
hy_0_snapshot[div(time,snap_step)] = (time, copy(hy_0).*globals.imp0)
end
end
In [5]:
anim = Animation()
for i = 1:num_snaps
p = plot(1:size, ez_snapshot[i][2], lab="Ez")
plot!(1:size, ez_0_snapshot[i][2], lab="Ez_0")
# plot!(1:size-1, hy_snapshot[i][2], lab="Hy*imp0")
time = ez_snapshot[i][1]
plot!(ann=[(0.8*size, 1.5, "time =$time")])
plot!([size/2, size/2], [-2, 2])
plot!([size/2+q_wavelength, size/2+q_wavelength], [-2, 2])
plot!(xlims=(1, size), ylims=(-2, 2))
frame(anim, p)
end
gif(anim, "./Task3/Quarter_Wavelength_Reflection.gif", fps=15)
Out[5]:
In [6]:
plot( (ez-ez_0)[1:div(size, 2)])
Out[6]:
In [ ]:
In [ ]: