% set up
s = sym('s','real');
theta(s) = sym('theta(s)');
rho(s) = sym('rho(s)');
z(s) = sym('z(s)');
k1(s) = sym('k1(s)');
k2(s) = sym('k2(s)');
k3(s) = sym('k3(s)');
eps_para = sym('eps_para','real');
eps_perp = sym('eps_perp','real');

% creates P_perp_k, projection onto plane perpendicular to k
k = [k1(s);k2;k3];
k_hat= k/sqrt(k1^2+k2^2+k3^2);
k_hat_otimes_k_hat = k_hat*k_hat.';
I = eye(3);
P_perp_k = I - k_hat_otimes_k_hat;
P_perp_k = simplify(P_perp_k);

% creates P_para_n of director
n = [cos(theta/2);sin(theta/2);0];
P_para_n = n*n.';
P_para_n = simplify(P_para_n);
P_perp_n = I - P_para_n;

% creates dielectric tensor
eps = eps_para*P_para_n+eps_perp*P_perp_n;

% creates Y matrix
Y = (k1^2 + k2^2 + k3^2)*P_perp_k*eps*P_perp_k;
Y = simplify(Y);

% calculates eigenvalues and simplifies them
E = simplify(eig(Y));

% identifies non-zero eigenvalue and removes from vector
for i=1:3,
    if E(i) == 0,
        E(i) = [];
        break
    end
end
clear i

% attempts to simplify E further
E = expand(E);
E = simplify(E, 'Steps', 1000);
E = factor(E);