% übernommen von
% philipp thaler
% 8.5.09

% solver für das Eigenwertproblem von 3d Dispersions-Relationen
function w = disp_solver(an, m, C, k)
% Richtungen der Nachbarn
ar =an(2:end,:);
arnorm = zeros(size(ar,1),1);
for n = 1:size(ar,1)
arnorm(n) = norm(ar(n,[1,2,3]));
ar(n,[1,2,3]) = ar(n,[1,2,3])./arnorm(n);
end
%Berechnung der Matrix laut Herleitung
for n = 1:size(ar,1)
% x-direction, sum(delta_u):
delta_uxx = (ar(:,1).*ar(:,1)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uxy = (ar(:,1).*ar(:,2)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uxz = (ar(:,1).*ar(:,3)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
% y-direction, sum(delta_u):
delta_uyx = (ar(:,2).*ar(:,1)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uyy = (ar(:,2).*ar(:,2)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uyz = (ar(:,2).*ar(:,3)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
% z-direction, sum(delta_u):
delta_uzx = (ar(:,3).*ar(:,1)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uzy = (ar(:,3).*ar(:,2)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
delta_uzz = (ar(:,3).*ar(:,3)).*(exp(i*sum(an(2:end,:).*repmat(k,size(ar,1),1),2))-1);
end
M = -real([sum(delta_uxx,1),sum(delta_uxy,1),sum(delta_uxz,1);sum(delta_uyx,1),sum(delta_uyy,1),sum(delta_uyz,1);sum(delta_uzx,1),sum(delta_uzy,1),sum(delta_uzz,1)]);
% eigs(M) = w^2 * m/C
%Berechnung von w laut hergeleiteter Formel
w = real(sqrt(eigs(M)*C/m));