% DISPRELELECTRON computes dispersion relations and root paths for waves on
% the electron time scale. 
% 
% [vph, UHPpole, lim, omega, f, k, v] = DISPRELELECTRON('inputfile',k,v)
%
% inputfile: m-file which defines the poles and residues of the distribution
% k:         wave number vector for which the solution is computed
% v:         velocity vector for which f is computed
%
% vph:     root paths in the complex phase velocity plane
% UHPpole: the poles of the distribution function in the upper half plane
% lim:     limiting solutions of vph for k=0
% omega:   complex angular frequency, Im{omega}>0 for stable solutions
% f:       distribution function

function [vph, UHPpole, lim, omega, f, k, v] = disprelelectron(inputfile, k, v)

% Defaults
if nargin<3
  v=[-3:0.005:3];
  if nargin<2
    k=[2.^(-6:0.1:8.3)];
  end
end

% Get the poles and residues
run(inputfile)

%%%% Limiting solution, k=0
P = polyelectron(UHPpole,UHPresid,0);            % Find disp.rel. polynomial

lim=roots(P);                                    % Solution to Eq.

sol_old=[1e6 -1e6 lim.'];                        % Initial value for the sort

omega=[]; vph=[];
for ii=1:length(k)
  P = polyelectron(UHPpole,UHPresid,k(ii)^2); % Find disp.rel. polynomial

% Issue a warning if the moduli of some coefficients are much larger than
% others
  if min(abs(P))/max(abs(P))<eps
    warning(['min(abs(P))/max(abs(P)) = ' ...
             num2str(min(abs(P))/max(abs(P)))])
  end

  solution = roots(P);               % This is where the equation is solved
 
  sol_sorted = refsort(sol_old, solution);%Sorts the roots into individual modes
  sol_old = sol_sorted;

  vph=[vph sol_sorted];
  omega = [omega k(ii)*sol_sorted];
end


% Distribution function
f=sum(ALLresid.'*ones(size(v))./ ...
       (ones(size(ALLpole)).'*v - ALLpole.'*ones(size(v))));
f=real(f);

