% POLYION finds the polynomial on the left-hand side of Eq. (8) of
% Gunell and Skiff Phys. Plasmas, vol. 8, 3550-3557 (2001). The
% right-hand side of that equation is 0. 
%
% P = POLYION(UHPpole, UHPresid, C)
%     finds the polynomial on the left-hand side of Eq. (8) of
%     Gunell and Skiff Phys. Plasmas, vol. 8, 3550-3557 (2001). The
%     right-hand side of that equation is 0.
%
% UHPpole:     the poles of the distribution function in the upper half plane
% UHPresidues: residues of those poles of the distribution function that
%              are in the upper half plane
% C:           a constant that should be equal to k^2. 
%
% P:           the polynomial that is the output of this function
%
%
% P = POLYION(UHPpole, UHPresid, C, ompd2, lde2)
%     also includes a component for which typical velocities are much
%     slower than the phase velocity, and which therefore can be
%     approximated by omega_pd^2/omega in Eq. (6) of the abovementioned
%     paper. This could for example be used for calculations of
%     dispersion relations for dust ion acoustic waves.
%
% ompd2:       square of the angular plasma frequency of the plasma
%              component approximated by omega_pd^2/omega.
% lde2:        square of the electron Debye length. If the electron Debye
%              length is used as reference for length, this should be set
%              to 1. Otherwise it should be expressed in units of the
%              reference length.

function P = polyion(UHPpole, UHPresid, C, ompd2, lde2)


if nargin<5
  lde2 = 1;
  if nargin<4
    ompd2 = 0;
  end
end

n=length(UHPpole);
fullpoly=poly(UHPpole);
fullpoly=conv(fullpoly,fullpoly); %% Double roots

if ompd2==0,
  P = (C+1/lde2)*fullpoly;
  for ii=1:n
    temp = deconv(fullpoly, poly(UHPpole(ii)));     % Polynomial division
    temp = deconv(temp    , poly(UHPpole(ii)));     % Polynomial division
    P    = P - 2i*pi*UHPresid(ii)*[0 0 temp];
  end
else
  P = (C+1/lde2)*[fullpoly 0 0] - ompd2*[0 0 fullpoly];
  for ii=1:n
    temp = deconv(fullpoly, poly(UHPpole(ii)));     % Polynomial division
    temp = deconv(temp    , poly(UHPpole(ii)));     % Polynomial division
    P    = P - 2i*pi*UHPresid(ii)*[0 0 temp 0 0];
  end
end
