Ran plasmonic couplers for longer, seems like efficiency increases.

Also started to correct power (Poynting vector) calculation.

paper_results.m

@@ -22,7 +22,8 @@
 % Create the specifications for the various design problems.
 specs = create_specs(dims, eps_lims, eps_uniform);
 
-for k = 1 : length(specs)
+% for k = 1 : length(specs)
+for k = 4 : 5
     eps{k} = solve(specs{k}, num_iters, 1e-6);
     simulate(specs{k}, eps{k}, [200 160]);
 end
@@ -94,7 +95,7 @@
 eps(end-1:end,  (dims(2)-w(2))/2:(dims(2)+w(2))/2) = 1; 
 eps(3:end-2, 3:end-2) = eps_uniform; % Fill the design area with uniform eps.
 
-spec{4} = setup(0.20, eps, eps_lims, [1 1]);
+spec{4} = setup(0.25, eps, eps_lims, [1 1]);
 
 
     %
@@ -107,4 +108,4 @@
 eps(end-1:end,  (dims(2)-w(2))/2:(dims(2)+w(2))/2) = -1.1; % Output wg.
 eps(3:end-2, 3:end-2) = eps_uniform; % Fill the design area with uniform eps.
 
-spec{5} = setup(0.20, eps, eps_lims, [1 1]);
+spec{5} = setup(0.25, eps, eps_lims, [1 1]);

simulate.m

@@ -27,7 +27,7 @@
 % Hard-coded parameters.
 t_pml = 10; % Thickness of PML.
 sigma_pml = 1 / spec.omega; % Strength of PML.
-exp_pml = 2.5; % Exponential spatial increase in pml strength.
+exp_pml = 3.5; % Exponential spatial increase in pml strength.
 
 
     % 
@@ -77,6 +77,7 @@
 
 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.
@@ -113,12 +114,29 @@
 out_pos = min([round(dims(1)-pad(2)/2), dims(1)-t_pml-1]); 
 
 % Project y onto x.
-proj = @(x, y) (dot(x(:), y(:)) / norm(x(:))^2) * x(:);
+proj = @(x, y) (dot(y(:), x(:)) / norm(x(:))^2) * x(:);
 
 % Calculate the power in the desired output mode.
-P_out = abs(dot(proj(spec.out.Hz, Hz(out_pos,pad(3)+1:end-pad(4))), ...
-                proj(spec.out.Ey * exp(i * 0.5 * spec.out.beta), ...
-                    Ey(out_pos,pad(3)+1:end-pad(4)))));
+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
 
 
     %