% DISPRELELECTRONNEW computes dispersion relations and root paths for waves on
% the electron time scale. 
%
% DISPRELELECTRONNEW does the same thing as DISPRELELECTRON. The
% difference is that DISPRELELECTRONNEW uses FMINSEARCH to find minima of
% |epsilon|^2, where epsilon is defined by Eq. (12) of Gunell and Skiff,
% Phys. Plasmas, vol. 9, 2585-2592 (2002). This means that we do not have
% to worry about the ill-posed problem of finding the roots to a
% polynomial where some coefficients are very small and some very
% large. Instead the difficulty is finding a good initial guess, and this
% may have to be fine tuned should this function be used for a particular
% problem. The solution in the k=0 limit is calculated in two
% ways. Neither of them is perfect. 
%
% [vph, UHPpole, lim,omega, f, k, v, lim2] = 
%        DISPRELELECTRONNEW('inputfile', k, v, opt, guess)
%
% 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
% opt:       structure containing the options for FMINSEARCH
% guess:     starting point for solution search
%
% 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
% lim2:    'limiting solutions' of vph computed for a small k>0

function [vph, UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew(inputfile, k, v, opt, guess)

% Defaults
if nargin<5
  guess=[];
  if nargin<4
    opt=optimset('MaxFunEvals',1000);
    if nargin<3
      v=[-3:0.005:3];
      if nargin<2
        k=[2.^(-6:0.1:8.3)];
      end
    end
  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.

% Start looking for solutions at sol_old
if isempty(guess)
  interpole=zeros(size(UHPpole));
  for ii=1:length(UHPpole);
    interpole(ii)=min(abs(UHPpole(ii)-UHPpole([1:ii-1 ii+1:end])));
  end
  sol_old = [(UHPpole+0.08*(1+1i)*abs(UHPpole)) ...
             (UHPpole-0.08*(1+1i)*abs(UHPpole))];
% $$$   sol_old = [(UHPpole+0.2*(1+1i)*interpole) ...
% $$$              (UHPpole-0.2*(1+1i)*interpole)];
else
  sol_old=guess;
end

omega=[]; vph=[];
for ii=length(k):-1:1    % We start at large k
  om1 = k(ii)*sol_old;
  for jj= 1:length(om1)
    x = fminsearch(@(om) epsilon2(om(1)+1i*om(2),k(ii),UHPpole,UHPresid),...
                   [real(om1(jj)) imag(om1(jj))],opt);
    om1(jj)=x(1)+1i*x(2);
  end
  disp(['countdown ' num2str(ii)])
  om1=om1(:);
  solution = om1/k(ii);
  sol_old = solution;

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

omega=fliplr(omega);
vph=fliplr(vph);


% Differently computed k=0 limit
om1 = 0.1*min(k)*sol_old;
for jj= 1:length(om1)
  x = fminsearch(@(om) epsilon2(om(1)+1i*om(2),0.1*min(k),UHPpole,UHPresid),...
                 [real(om1(jj)) imag(om1(jj))],opt);
  om1(jj)=x(1)+1i*x(2);
end
lim2 = om1 / (0.1*min(k));

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


% -------------------------------------------------- %
function ep2 = epsilon2(omega, k, UHPpole, UHPresid)
  ep=1-2i*pi*sum(UHPresid./(omega-k*UHPpole).^2);
  ep2=ep*ep';
