ob1_wgmode.m

function [beta, Hz, Ey] = ob1_wgmode(omega, eps, order)
% Waveguide mode solver.
% 
% Casts the problem in terms of Hz components, solves for eigenvalues and 
% eigenvectors, and then selects the correct order of eigenvalue.
% Assumes that the mode is propagating in the x-direction.
% 
% If ORDER = 0, then asks user to select one of the plotted modes.
% 
% If user has selected an evanescent mode, then gives a warning.

N = length(eps); 


    % 
    % Setup the eigenvalue problem.
    %

% Get the (smoothed) x and y components of epsilon.
[eps_x, eps_y] = ob1_interp_eps(eps); 

% Build the derivative matrix.
Dy = spdiags([-ones(N,1), ones(N,1)], [-1 0], N, N);
Dy(1,N) = -1;

% Build the matrix for the eigenvalue problem.
my_diag = @(x) spdiags(x(:), 0, N, N);
A = my_diag(eps_y) * (-Dy' * my_diag(eps_x.^-1) * Dy + omega^2 * speye(N));  

% % Remove the wrap-around elements (we don't want periodic boundaries).
% A(N,1) = 0;
% A(1,N) = 0;


   %
   % Solve the eigenvalue problem.
   %

[V, D] = eig(full(A)); % Just use the dense eigenvalue solver.

% Sort modes based on descending beta^2.
[beta2, ind] = sort(real(diag(D)), 'descend');
beta = sqrt(beta2);

% Sort the mode profiles.
Hz = V(:,ind);

% Orient modes so that the largest magnitude element is always positive.
Hz = Hz * diag((max(Hz) > max(-Hz)) - (max(Hz) <= max(-Hz)));

% Calculate Ey.
Ey = ((1 / omega) * my_diag(eps_y(:).^-1) * Hz) * my_diag(beta);

% Normalize power (to unity).
power = diag(Ey' * Hz);
Hz = Hz * diag(power.^-0.5);
Ey = Ey * diag(power.^-0.5);


    %
    % Choose the desired mode.
    %

if (order == 0) % Cycle through modes one at a time and let user choose.
    fprintf('Enter "" (blank) to view next mode, ');
    fprintf('or the number of the mode to be selected.\n');

    for k = 1 : N
        my_mode_plot(Hz(:,k), Ey(:,k)); % Let user see the mode.
        prompt = sprintf('mode %d: ', k); % Tell the user which mode this is.
        choice = input(prompt, 's'); % Ask the user if they want it.

        if strcmp(choice, '') % User wants to see next mode.
            % Delete previous line so we can write on-top of it.
            fprintf(repmat('\b', 1, length(prompt) + length(choice) + 1));
        else % User has chosen a mode
            order = str2num(choice);
            beta = beta(order);
            Hz = Hz(:,order);
            Ey = Ey(:,order);
            break
        end
    end

else % User has already specified what mode they want, so just choose it.
    beta = beta(order);
    Hz = Hz(:,order);
    Ey = Ey(:,order);
end

if ~isreal(beta) % Not a propagating mode, warn user!
    warning('Selected mode is not propagating (non-real wave-vector)!');
end

    %
    % Print out results.
    %

err = norm(A * Hz - beta^2 * Hz);
power = Ey' * Hz;

fprintf('Mode number: %d\n', order);
fprintf('Beta (wave-vector): %1.3f\n', beta);
fprintf('Error in eigenvalue equation: %1.1e\n', err);
fprintf('Power: %1.3f\n\n', power);

% Plot result.
my_mode_plot(Hz, Ey);



function my_mode_plot(Hz, Ey)
% Visualize both components of the waveguide mode.
subplot 121; my_plot(Hz(:)); title('Hz');
subplot 122; my_plot(Ey(:)); title('Ey');

function my_plot(z)
% Plotting function.
area(1:length(z), z)
axis([1 length(z) -1.5*max(abs(z)) 1.5*max(abs(z))])