%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the propagation constant for a general
% material. The direction of propagation is given by
% theta and phi, and the frequency is fixed.
%
% Daniel Sjoberg, 2012-09-27
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n, E, H] = propconst2(theta, phi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Propagation parameters and constants
%
f = 10e9; % Frequency in Hz
c0 = 299792458; % Speed of light in vacuum
mu0 = 4*pi*1e-7; % Permeability of vacuum
eps0 = 1/c0^2/mu0; % Permittivity of vacuum
eta0 = sqrt(mu0/eps0); % Wave impedance of vacuum
omega = 2*pi*f; % Angular frequency
k0 = omega/c0; % Wave number in vacuum
kt = [0; 0]; % Transverse wave vector = 0,
% propagation is purely in the
% z direction in the material
% coordinate system
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Material parameters
%
epsilon = eps0*[1, 0, 0
0, 2, 0
0, 0, 3];
xi = 1/c0*[1, 0, 0
0, 1, 0
0, 0, 1]*0;
zeta = 1/c0*[1, 0, 0
0, 1, 0
0, 0, 1]*0;
mu = mu0*[1, 0, 0
0, 1, 0
0, 0, 1];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Rotate the material matrices so that the propagation
% direction indicated by theta and phi corresponds to the
% z direction in the material coordinate system
%
Rphi = [cos(phi), -sin(phi), 0
sin(phi), cos(phi), 0
0, 0, 1];
Rtheta = [cos(theta), 0, sin(theta)
0, 1, 0
-sin(theta), 0, cos(theta)];
R = Rtheta*Rphi;
epsilon = R*epsilon*transpose(R);
xi = R*xi*transpose(R);
zeta = R*zeta*transpose(R);
mu = R*mu*transpose(R);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Splitting of the material matrices into submatrices
% and vectors.
%
epsilontt = epsilon(1:2, 1:2);
epsilont = epsilon(1:2, 3);
epsilonz = epsilon(3, 1:2);
epsilonzz = epsilon(3, 3);
xitt = xi(1:2, 1:2);
xit = xi(1:2, 3);
xiz = xi(3, 1:2);
xizz = xi(3, 3);
zetatt = zeta(1:2, 1:2);
zetat = zeta(1:2, 3);
zetaz = zeta(3, 1:2);
zetazz = zeta(3, 3);
mutt = mu(1:2, 1:2);
mut = mu(1:2, 3);
muz = mu(3, 1:2);
muzz = mu(3, 3);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the A matrix
%
zcross = [0, -1
1, 0];
Mzz = [epsilonzz, xizz
zetazz, muzz];
ktprim = zcross*kt;
Zv = [0
0];
A1 = [Zv, ktprim/omega
-ktprim/omega, Zv] - [epsilont, xit
zetat, mut];
A2 = [transpose(Zv), -transpose(ktprim/omega)
transpose(ktprim/omega), transpose(Zv)] - [epsilonz, xiz
zetaz, muz];
A = A1*inv(Mzz)*A2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Finally, put the fundamental matrix together
%
Mtt = [epsilontt, xitt
zetatt, mutt];
% Using [Et,Htxz] as eigenvectors
zleft = [0, 0, 0, 1
0, 0, -1, 0
1, 0, 0, 0
0, 1, 0, 0];
zright = [1, 0, 0, 0
0, 1, 0, 0
0, 0, 0, -1
0, 0, 1, 0];
W = zleft*(Mtt - A)*zright;
% Uncomment the following lines to use [Et,Ht] as eigenvectors
% zcrossJ = [0, 0, 0, 1
% 0, 0, -1, 0
% 0, -1, 0, 0
% 1, 0, 0, 0];
% W = zcrossJ*(Mtt - A);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the eigenvalues and eigenvectors
%
[v, ntmp] = eig(W);
n_over_c0 = diag(ntmp); % Pick out the diagonal elements
% First sort the refractive indices in descending order according
% to the real part. The first two elements will correspond to the
% waves propagating in the positive z direction.
[tmp, idx] = sort(real(n_over_c0), 1, 'descend');
E = v(1:2,idx(1:2));
H = [-v(4,idx(1:2)) % This applies if [Et,Htxz] are eigenvectors
v(3,idx(1:2))]; % If [Et,Ht] are eigenvectors, H=v(3:4,:);
n = n_over_c0(idx(1:2))*c0;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%