From 16ddf8c4f0daa631ee5c5447d47ebe0f23112b21 Mon Sep 17 00:00:00 2001
From: Jesse Lu <mr.jesselu@gmail.com>
Date: Mon, 24 Oct 2011 16:43:45 -0700
Subject: [PATCH] Draft of power calculation works.

---
 ob1_calc_power.m | 104 +++++++++++++++++++++++++++++++++--------------
 ob1_fdfd.m       |  48 +++++++++++++++-------
 ob1_pad_eps.m    |   5 +++
 simulate.m       |  41 ++++++++++++-------
 4 files changed, 138 insertions(+), 60 deletions(-)
 create mode 100644 ob1_pad_eps.m

diff --git a/ob1_calc_power.m b/ob1_calc_power.m
index 834cf47..f6f1fcc 100644
--- a/ob1_calc_power.m
+++ b/ob1_calc_power.m
@@ -1,38 +1,80 @@
-function [P_out] = ob1_calc_power(spec, Ex, Ey, Hz, t_pml, pad)
+function [Pout, err] = ob1_calc_power(E, H, mode)
+% Description
+%     Calculate power in an output mode. Assumes that all the power in the
+%     mode is propagating to the right (no reflections).
 
-dims = size(Ex);
+dims = size(E);
+
+    
+    % 
+    % Obtain magnitude of mode in both E and H fields.
     %
-    % Calculate power output to desired mode.
+
+E_mag = calc_mag(mode.Ey, E);
+H_mag = calc_mag(mode.Hz, H);
+
+fields = [E_mag(:), H_mag(:)];
+
+    
+    % 
+    % Estimate the amplitude, wave-vector and phase of the 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(:);
+phase = unwrap(angle(fields));
+x = [[1:dims(1)]-0.5; [1:dims(1)]]'; % Accounts for offset in Yee grid.
+
+A = mean(abs(fields), 1); 
+[beta(1), phi(1)] = my_linear_regression(x(:,1), phase(:,1));
+[beta(2), phi(2)] = my_linear_regression(x(:,2), phase(:,2));
+
+
+    %
+    % Local optimization to tune fit parameters.
+    %
+
+p = mean([A; beta; phi], 2);
+fun = @(p) p(1) * exp(i * (p(2) * x + p(3)));
+err_fun = @(p) (1/norm(fields, 'fro')) * ...
+        norm(fields - fun(p));
+options = optimset('MaxFunEvals', 1e3);
+p = fminsearch(err_fun, p, options); % Optimize.
+
+
+    %
+    % Calculate output power and error.
+    %
 
-% 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)))));
+amplitude = p(1);
+Pout = amplitude^2;
 
-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));
+err = err_fun(p);
+err_limit = 1e-3;
+if (err > err_limit) % If error is somewhat large, tell user.
+    warning('Error in fit exceeds threshold (%e > %e).', err, err_limit);
 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
-  
+
+% Plot.
+plot(real(fields), 'o-');
+hold on
+plot(real(fun(p)), '.-');
+hold off
+
+
+function [mag] = calc_mag (ref, field)
+% Project y onto x.
+my_dot = @(x, y) (dot(y(:), x(:)) / norm(x(:))^2);
+
+for k = 1 : size(field, 1)
+    mag(k) = my_dot(ref, field(k,:));
+end
+
+
+
+function [a, b] = my_linear_regression(x, y)
+% Fit data to simple model, y = a * x + b
+y_ctr = y - mean(y);
+x_ctr = x - mean(x);
+a = (x_ctr' * y_ctr) / norm(x_ctr)^2;
+b = -(a * mean(x) - mean(y));
+
+
diff --git a/ob1_fdfd.m b/ob1_fdfd.m
index f84f72a..c92bfa6 100644
--- a/ob1_fdfd.m
+++ b/ob1_fdfd.m
@@ -1,15 +1,31 @@
-function [Ex, Ey, Hz] = ob1_fdfd(spec, eps, t_pml, pad)
+function [Ex, Ey, Hz] = ob1_fdfd(omega, eps, in)
 % 
 % Description
-%     Solve an FDFD (finite-difference, frequency-domain) problem.
+%     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.
+
 
-sigma_pml = 1 / spec.omega; % Strength of PML.
+% Hard-coded parameters for the pml.
+t_pml = 10; % Thickness of pml.
+sigma_pml = 1 / 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.
 
+% Expand eps to include room for the pml padding
+eps = ob1_pad_eps(eps, t_pml * [1 1 1 1]);
+dims = size(eps); % New size of the structure.
+[eps_x, eps_y] = ob1_interp_eps(eps); % Obtain x- and y- components of eps.
+
+    
     %
-    % Build the simulation matrix.
+    % Form the system matrix which will be solved.
     %
 
 % Shortcut to form a derivative matrix.
@@ -29,7 +45,7 @@
 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));
+A = Ecurl * inv_eps * Hcurl - omega^2 * speye(prod(dims));
 
     
     %
@@ -37,19 +53,17 @@
     %
 
 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;
+in_pos = 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.
-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(in_pos+1, pad(3)+1:end-pad(4)) = in.Hz;
+b_pad = [floor((dims(2) - length(in.Hz))/2), ...
+        ceil((dims(2) - length(in.Hz))/2)];
+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.
 
-% 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.
 
 
@@ -59,10 +73,14 @@
 
 Hz = A \ b; % This should be using sparse matrix factorization. 
 
-E = 1/spec.omega * inv_eps * Hcurl * Hz; % Obtain E-fields.
+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.
+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);
diff --git a/ob1_pad_eps.m b/ob1_pad_eps.m
new file mode 100644
index 0000000..2e39c44
--- /dev/null
+++ b/ob1_pad_eps.m
@@ -0,0 +1,5 @@
+function [eps] = ob1_pad_eps(eps, pad)
+
+% Expand eps to the full simulation size.
+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)));
diff --git a/simulate.m b/simulate.m
index cfba8d0..c92d791 100644
--- a/simulate.m
+++ b/simulate.m
@@ -24,8 +24,6 @@
 %         an unbroken version of the input waveguide, this is the true
 %         (accurate) amount of power excited in the input mode.
 
-% Hard-coded parameters.
-t_pml = 10; % Thickness of PML.
 
     
     % 
@@ -38,28 +36,43 @@
         floor((dims(2) - size(eps, 2))/2), ...
         ceil((dims(2) - size(eps, 2))/2)];
 
-% Expand eps to the full simulation size.
-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 = ob1_pad_eps(eps, pad);
 
 
-[Ex, Ey, Hz] = ob1_fdfd(spec, eps, t_pml, pad); % Solve.
+    %
+    % Determine input power by solving a structure with only an input
+    % waveguide.
+    %
+
+[Ex, Ey, Hz] = ob1_fdfd(spec.omega, repmat(eps(1,:), size(eps,1), 1), spec.in);
+
+x_out = round(dims(1)-pad(2)/2+1) : dims(1);
+y_out = pad(3)+1 : dims(2)-pad(4);
+P_in = ob1_calc_power(Ey(x_out,y_out), Hz(x_out,y_out), spec.in);
+
+
+    %
+    % Calculate output power for the structure (and actual output mode).
+    %
+
+[Ex, Ey, Hz] = ob1_fdfd(spec.omega, eps, spec.in);
+P_out = ob1_calc_power(Ey(x_out,y_out), Hz(x_out,y_out), spec.out);
+
 
-[P_out] = ob1_calc_power(spec, Ex, Ey, Hz, t_pml, pad);
     %
     % Print and plot results.
     %
 
-fprintf('Output power in desired mode (input power approx. 1.0) : %1.3f\n', ...
-    P_out);
+fprintf('P_in: %1.3e \nP_out: %1.3e \neff: %1.3f%%\n', ...
+    P_in, P_out, 100 * P_out / P_in);
 
 ob1_plot(dims, {'\epsilon', eps}, {'|Hz|', abs(Hz)}, {'Re(Hz)', real(Hz)});
 
-% % The following commands may be used (uncommented) in order to plot more
-% % field information.
-% % figure(1); 
-% ob1_plot(dims, {'\epsilon', eps}, {'|Hz|', abs(Hz)}, {'Re(Hz)', real(Hz)});
-% 
+% The following commands may be used (uncommented) in order to plot more
+% field information.
+% figure(1); 
+ob1_plot(dims, {'\epsilon', eps}, {'|Hz|', abs(Hz)}, {'Re(Hz)', real(Hz)});
+
 % % Plot all fields.
 % figure(2); 
 % ob1_plot(dims, ...