% S3Dofk Spectral density as a function of wave number k and angular
% frequency omega for three-dimensional isotropic distributions. 
%
% S = S3Dofk(k, omega, ALLpole, ALLresid, UHPpole, UHPresid, ...
%            n_vec, nb_vec, mass_vec, me)
%
% ALLpole, ALLresid, UHPpole, UHPresid, n_vec, nb_vec describe the
% projection of the distribution function on an axis parallel to k. 
%
% The mass_vec vector specifies the mass for the particle species
% corresponding to each element in n_vec and nb_vec.
%
% This file is meant to be integrated according to Eq. (14) of 
% Gunell and Skiff, Phys. Plasmas, vol. 9, 2585-2592 (2002).
% To obtain S as defined in Eq. (13) of that paper use S1Dofk instead.

function S = S3Dofk(k, omega, ALLpole, ALLresid, UHPpole, UHPresid, ...
                   n_vec, nb_vec, mass_vec, me)

% Phase velocity
v=omega./k;

% Initialise
epsilon=zeros(size(k));
S=epsilon;

% Set up masses vector for use below
masses=[];
for ii=1:length(mass_vec)
  masses = [masses mass_vec(ii)*ones(1,2*(n_vec(ii)+nb_vec(ii)))];
end

% Compute epsilon and the modulus of it squared
for jj=1:length(UHPpole)
  epsilon=epsilon + UHPresid(jj)./(omega-k*UHPpole(jj)).^2;
end
epsilon=1-2i*pi*epsilon;
epsilon2=abs(epsilon).^2;

% Form S by summation over all poles
for ii=1:length(ALLpole)
  S = S + masses(ii)*ALLresid(ii)./(v-ALLpole(ii));
end

% And finally
S = real(pi*S./(k.*epsilon2))/me;

