diff --git a/ob1_calc_power.m b/ob1_calc_power.m new file mode 100644 index 0000000..834cf47 --- /dev/null +++ b/ob1_calc_power.m @@ -0,0 +1,38 @@ +function [P_out] = ob1_calc_power(spec, Ex, Ey, Hz, t_pml, pad) + +dims = size(Ex); + % + % Calculate power output to desired mode. + % + +% Location for power calculation. +out_pos = min([round(dims(1)-pad(2)/2), dims(1)-t_pml-1]); + +% Project y onto x. +proj = @(x, y) (dot(y(:), x(:)) / norm(x(:))^2) * x(:); + +% Calculate the power in the desired output mode. +calcP = @(loc) 0.5 * real(... + dot(proj(spec.out.Hz, Hz(loc,pad(3)+1:end-pad(4))), ... + proj(spec.out.Ey * exp(0.0 * i * spec.out.beta), ... + Ey(loc,pad(3)+1:end-pad(4))))); + +out_pos = round(dims(1) - pad(2)) : dims(1) - t_pml - 1; +for k = 1 : length(out_pos) + P_out(k) = calcP(out_pos(k)); +end +plot(P_out, '.-') +P_out +pause +P_out = mean(P_out) + +% Calculate power leaving a box. +Pbox = @(x,y) dot(Hz(x,y), Ey(x,y)); +box_pad = t_pml + 5; +box = [box_pad, dims(1)-box_pad, box_pad, dims(2)-box_pad]; +% bottom = 0.5 * real(Pbox(box(1):box(2),box(3))) +% top = 0.5 * real(Pbox(box(1):box(2),box(4))) +% left = 0.5 * real(Pbox(box(1),box(3):box(4))) +right = 0.5 * real(Pbox(box(2),box(3):box(4))) +% Pbox_total = bottom + top + left + right + diff --git a/ob1_fdfd.m b/ob1_fdfd.m new file mode 100644 index 0000000..f84f72a --- /dev/null +++ b/ob1_fdfd.m @@ -0,0 +1,68 @@ +function [Ex, Ey, Hz] = ob1_fdfd(spec, eps, t_pml, pad) +% +% Description +% Solve an FDFD (finite-difference, frequency-domain) problem. + +sigma_pml = 1 / spec.omega; % Strength of PML. +exp_pml = 3.5; % Exponential spatial increase in pml strength. +dims = size(eps); +[eps_x, eps_y] = ob1_interp_eps(eps); % Get x and y components of eps. + + % + % Build the simulation matrix. + % + +% 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. +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))]; + +% 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 - spec.omega^2 * speye(prod(dims)); + + + % + % Determine the input excitation. + % + +b = zeros(dims); % Input excitation, equivalent to magnetic current source. +% in_pos = max([t_pml+1, round(pad(1)/2)]); % Location of input excitation. +in_pos = 5; + +% For one-way excitation in the forward (to the right) direction, +% we simple cancel out excitation in the backward (left) direction. +b(in_pos+1, pad(3)+1:end-pad(4)) = spec.in.Hz; +b(in_pos, pad(3)+1:end-pad(4)) = -spec.in.Hz * exp(i * spec.in.beta); + +b = b ./ eps_y; % Convert from field to current source. + +% Normalization factor so that the input power is unity. +b = -i * 2 * spec.in.beta / (1 - exp(i * 2 * spec.in.beta)) * b; + +b = b(:); % Vectorize. + + + % + % Solve. + % + +Hz = A \ b; % This should be using sparse matrix factorization. + +E = 1/spec.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); + diff --git a/simulate.m b/simulate.m index f47702d..cfba8d0 100644 --- a/simulate.m +++ b/simulate.m @@ -26,8 +26,6 @@ % Hard-coded parameters. t_pml = 10; % Thickness of PML. -sigma_pml = 1 / spec.omega; % Strength of PML. -exp_pml = 3.5; % Exponential spatial increase in pml strength. % @@ -44,101 +42,10 @@ eps = cat(1, repmat(eps(1,:), pad(1), 1), eps, repmat(eps(end,:), pad(2), 1)); eps = cat(2, repmat(eps(:,1), 1, pad(3)), eps, repmat(eps(:,end), 1, pad(4))); -[eps_x, eps_y] = ob1_interp_eps(eps); % Get x and y components of eps. +[Ex, Ey, Hz] = ob1_fdfd(spec, eps, t_pml, pad); % Solve. - % - % Build the simulation matrix. - % - -% 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. -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))]; - -% 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 - spec.omega^2 * speye(prod(dims)); - - - % - % Determine the input excitation. - % - -b = zeros(dims); % Input excitation, equivalent to magnetic current source. -in_pos = max([t_pml+1, round(pad(1)/2)]); % Location of input excitation. -in_pos = 5; - -% For one-way excitation in the forward (to the right) direction, -% we simple cancel out excitation in the backward (left) direction. -b(in_pos+1, pad(3)+1:end-pad(4)) = spec.in.Hz; -b(in_pos, pad(3)+1:end-pad(4)) = -spec.in.Hz * exp(i * spec.in.beta); - -b = b ./ eps_y; % Convert from field to current source. - -% Normalization factor so that the input power is unity. -b = -i * 2 * spec.in.beta / (1 - exp(i * 2 * spec.in.beta)) * b; - -b = b(:); % Vectorize. - - - % - % Solve. - % - -Hz = A \ b; % This should be using sparse matrix factorization. - -E = 1/spec.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); - - - % - % Calculate power output to desired mode. - % - -% Location for power calculation. -out_pos = min([round(dims(1)-pad(2)/2), dims(1)-t_pml-1]); - -% Project y onto x. -proj = @(x, y) (dot(y(:), x(:)) / norm(x(:))^2) * x(:); - -% Calculate the power in the desired output mode. -calcP = @(loc) 0.5 * real(... - dot(proj(spec.out.Ey, Ey(out_pos,pad(3)+1:end-pad(4))), ... - proj(spec.out.Hz, Hz(out_pos,pad(3)+1:end-pad(4))))); - -out_pos = round(dims(1) - pad(2)) : dims(1) - t_pml - 1; -for k = 1 : length(out_pos) - P_out(k) = calcP(out_pos(k)); -end -plot(P_out, '.-') -P_out = mean(P_out) - -% Calculate power leaving a box. -Pbox = @(x,y) dot(Ey(x,y), Hz(x,y)); -box_pad = t_pml + 5; -box = [box_pad, dims(1)-box_pad, box_pad, dims(2)-box_pad]; -bottom = 0.5 * real(Pbox(box(1):box(2),box(3))) -top = 0.5 * real(Pbox(box(1):box(2),box(4))) -left = 0.5 * real(Pbox(box(1),box(3):box(4))) -right = 0.5 * real(Pbox(box(2),box(3):box(4))) -Pbox_total = bottom + top + left + right - - +[P_out] = ob1_calc_power(spec, Ex, Ey, Hz, t_pml, pad); % % Print and plot results. %