% aa16illustration input file for disprelion.
% This file is used to illustrate how simple pole expansion distributions
% can be constructed.
%
% HG 2016-09-08


ALLpole =[];
ALLresid=[];   

n_vec    =[2 0];                        % # of terms for Maxwell expansion
nb_vec   =[0 3];                        % # of terms for Butterworth
op2      =[1 1];                        % density
vtbut_vec=[1.5 1.5];                    % Cut velocity in cs units
vdbut_vec=[0.3 0.3];                    % Butterworth center velocity
vt       =[0.55 0.55];                  % thermal velocity, cs=1
vd       =[2.8 2.8];                    % drift velocity, cs=1


%op2=(op2/sum(op2)); % We don't normalise quite in the normal way now,
                     % because we want to plot M(v) and T(v) separately.
%%%% Normalise
vtbut_vec=vtbut_vec./vt;
vdbut_vec=vdbut_vec./vt;
for jj=1:length(n_vec)
  n= n_vec(jj);  
  nb=nb_vec(jj);
  vtbut=vtbut_vec(jj);
  vdbut=vdbut_vec(jj);

  pol=0*(1:2*n+1);
  pol(2*n+1)=1;
  for ii=1:n                               % Poles for Maxwellian
    pol(2*(n-ii)+1)=1.0/(2^i*prod(1:ii));
  end
  sol=roots(pol);
  res=pole2residue(sol);                   % Find residues

  z0=butterworthpoles(nb,vtbut); 	   % Butterworth poles
  z0=z0+vdbut;

  allpole =[sol;z0];			   % poles from expansion and mask
  allresid=pole2residue(allpole);          % Find residues
  [uhppole,uhpresid]=uhp(allpole,allresid);% Keep poles in upper half plane

  allresid=allresid/(2i*pi*sum(uhpresid)); % normalise zeroth moment
  for kk=1:length(allpole)                 % return to variables with dimension
    ALLpole =[ALLpole allpole(kk)*vt(jj)+vd(jj)];
    ALLresid=[ALLresid op2(jj)*allresid(kk)];
  end
  
end

[UHPpole,UHPresid]=uhp(ALLpole,ALLresid); % Keep poles in upper half plane
