% [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