% [w,U,x] = waveguide(k, n, M, L, h, a1, b1)
%
% Usage: compute the n smallest-frequency (w) waveguide modes
% of b div(a grad u) = -w^2 u, where u = exp(iky) u_k(x)
% and b,a = b1,a1 in a region of width h around x=0,
% and = 1 otherwise.
%
% Use a computational cell of M points and total width L in the
% x direction, with Dirichlet boundary conditions. Note that
% the boundary conditions are of negligible importance for the
% guided modes if L is big enough, since the guided modes are
% exponentially localized. The guided modes are the ones for
% which w < k.
%
% Return a row vector of n eigenvalues w, and a matrix U whose
% n columns are the eigenvectors, for the column vector x of
% coordinates.
%
% k can be a vector, in which case w is a length(k) by n matrix
% of eigenvalues at each k, and U are the eigenvectors for k(end).
function [w,U,x] = waveguide(k, n, M, L, h, a1, b1)
dx = L / (M+1);
x = [1:M]'*dx - L/2;
a = 1 + (a1 - 1) * (abs(x) < h/2);
b = 1 + (b1 - 1) * (abs(x) < h/2);
x2 = [0:M]'*dx + dx/2 - L/2;
a2 = 1 + (a1 - 1) * (abs(x2) < h/2);
D = sdiff1(M);
A = spdiags(-b,0,M,M) * D' * spdiags(a2,0,M+1,M+1) * D;
U = zeros(length(x), n);
w = zeros(length(k), n);
for i = 1:length(k)
Ak = A - k(i)^2 * spdiags(a.*b, 0, M, M);
[U,S] = eigs(Ak, n, 'sm');
[lam,ilam] = sort(-diag(S)); % make sure frequencies are in order
w(i,:) = sqrt(lam);
U = U(:,ilam);
% eigs normalizes vectors to length 1, but picks a random sign
% ... it is convenient to pick a deterministic sign based
% on average U in first half (1:M/2).
for j = 1:n
U(:,j) = U(:,j) * sign(mean(U(1:round(M/2),j)));
end
end