%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the propagation constant for a general
% material. The direction of propagation is z, and
% the frequency is fixed.
%
% Daniel Sjoberg, 2012-09-10
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Propagation parameters and constants
%
f = 10e9; % Frequency in Hz
theta = 0*pi/180; % The transverse wave vector is
ktdir = [1 % later represented as
0]; % kt=k0*sin(theta)*ktdir
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 = k0*sin(theta)*ktdir; % Transverse wave vector
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Material parameters
%
epsilon = eps0*[1, 0, 0
0, 1, 0
0, 0, 1]*(2-0.1j);
xi = 1/c0*[0, 0, 0
0, 0, 0
0, 0, 0];
zeta = 1/c0*[0, 0, 0
0, 0, 0
0, 0, 0];
mu = mu0*[1, 0, 0
0, 1, 0
0, 0, 1];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 and
% display the results
%
[v, n_over_c0] = eig(W);
E = v(1:2,:);
H = [-v(4,:) % This applies if [Et,Htxz] are eigenvectors
v(3,:)]; % If [Et,Ht] are eigenvectors, H=v(3:4,:);
n = n_over_c0*c0;
disp('Eigenvalues (refractive index):')
disp(n)
disp('Eigenvectors (polarization):')
disp(v)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%