-
Notifications
You must be signed in to change notification settings - Fork 18
/
ob1_fdfd.m
105 lines (83 loc) · 3.34 KB
/
ob1_fdfd.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
function [Ex, Ey, Hz] = ob1_fdfd(omega, eps, in, bc)
%
% Description
% Solve a FDFD (finite-difference, frequency-domain) problem, using the
% input mode as the source term.
%
% OB1_FDFD actually expands the structure in order to add absorbing
% boundary layers (pml), simulates the structure by sourcing it from
% within the pml, and then truncates the solution so the user does not
% see the absorbing pml region.
%
% Note that a time dependence of exp(-i * omega * t) is assumed for all
% fields.
% Hard-coded parameters for the pml.
t_pml = 20; % Thickness of pml.
sigma_pml = 1 / omega; % Strength of pml.
exp_pml = 2.5; % Exponential spatial increase in pml strength.
% Expand eps to include room for the pml padding
if strcmp(bc, 'pml')
eps = ob1_pad_eps(eps, t_pml * [1 1 1 1]);
elseif strcmp(bc, 'per')
eps = ob1_pad_eps(eps, t_pml * [1 1 0 0]);
end
dims = size(eps); % New size of the structure.
[eps_x, eps_y] = ob1_interp_eps(eps); % Obtain x- and y- components of eps.
%
% Form the system matrix which will be solved.
%
% Shortcut to form a derivative matrix.
S = @(sx, sy) ob1_shift_matrix(dims, -[sx sy]);
% Helper function to create stretched-coordinate PML absorbing layers.
scx = @(sx, sy) ob1_stretched_coords(dims, [1 dims(1)+0.5], [sx, sy], ...
'x', t_pml, sigma_pml, exp_pml);
scy = @(sx, sy) ob1_stretched_coords(dims, [1 dims(2)+0.5], [sx, sy], ...
'y', t_pml, sigma_pml, exp_pml);
% Define the curl operators as applied to E and H, respectively.
if strcmp(bc, 'pml')
Ecurl = [scy(.5,.5)*-(S(0,1)-S(0,0)), scx(.5,.5)*(S(1,0)-S(0,0))];
Hcurl = [scy(.5,0)*(S(0,0)-S(0,-1)); scx(0,.5)*-(S(0,0)-S(-1,0))];
elseif strcmp(bc, 'per')
Ecurl = [-(S(0,1)-S(0,0)), scx(.5,.5)*(S(1,0)-S(0,0))];
Hcurl = [(S(0,0)-S(0,-1)); scx(0,.5)*-(S(0,0)-S(-1,0))];
end
% Diagonal matrix for 1/epsilon.
inv_eps = spdiags([eps_x(:).^-1; eps_y(:).^-1], 0, 2*prod(dims), 2*prod(dims));
% This is the matrix that we will solve.
A = Ecurl * inv_eps * Hcurl - omega^2 * speye(prod(dims));
%
% Determine the input excitation.
%
b = zeros(dims); % Input excitation, equivalent to magnetic current source.
in_pos = t_pml+2; % Cannot be 1, because eps interpolation wreaks havoc at border.
% For one-way excitation in the forward (to the right) direction,
% we simple cancel out excitation in the backward (left) direction.
if strcmp(bc, 'pml')
b_pad = [floor((dims(2) - length(in.Hz))/2), ...
ceil((dims(2) - length(in.Hz))/2)];
elseif strcmp(bc, 'per')
b_pad = [0 0];
end
b(in_pos+1, b_pad(1)+1:end-b_pad(2)) = in.Hz;
b(in_pos, b_pad(1)+1:end-b_pad(2)) = -in.Hz * exp(i * in.beta);
b = b ./ eps_y; % Convert from field to current source.
b = b(:); % Vectorize.
%
% Solve.
%
Hz = A \ b; % This should be using sparse matrix factorization.
E = (i/omega) * inv_eps * Hcurl * Hz; % Obtain E-fields.
% Reshape and extract all three fields.
Ex = reshape(E(1:prod(dims)), dims);
Ey = reshape(E(prod(dims)+1:end), dims);
Hz = reshape(Hz, dims);
% Cut off the pml parts.
if strcmp(bc, 'pml')
Ex = Ex(t_pml+1:end-t_pml, t_pml+1:end-t_pml);
Ey = Ey(t_pml+1:end-t_pml, t_pml+1:end-t_pml);
Hz = Hz(t_pml+1:end-t_pml, t_pml+1:end-t_pml);
elseif strcmp(bc, 'per')
Ex = Ex(t_pml+1:end-t_pml, :);
Ey = Ey(t_pml+1:end-t_pml, :);
Hz = Hz(t_pml+1:end-t_pml, :);
end