% fluct1 distribution with both protons and electrons for fluctuation
% calculations. 
% This is meant to reproduce the distribution on the last line of Table I
% in the article by Gunell and Skiff, Phys. Plasmas, vol. 9, 2585-2592
% (2002).

ALLpole =[];
ALLresid=[];

me=9.10938291e-31;
mp=1.672621777e-27;

mass_vec =[mp mp me];                         % particle masses in kg

n_vec    =[5 3 3];                            % terms Maxwell expansion
nb_vec   =[0 0 0];                            % Butterworth 
vtbut_vec=[0 0 0]/sqrt(mp/me);                % Cutoff velocity in vte units
vdbut_vec=[0 0 0];                            % Butterworth center velocity
op2      =[0.4*me/mp 0.6*me/mp 1];            % plasma frequencies
vt       =[0.1 0.4  sqrt(mp/me)]/sqrt(mp/me); % thermal velocity
vd       =[0.0 0.0  0.0]/sqrt(mp/me);         % drift velocity vte==1


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^ii*prod(1:ii));
  end
  sol=roots(pol);
  res=pole2residue(sol);

  z0=butterworthpoles(nb,vtbut);
  z0=z0+vdbut;

  allpole =[sol;z0];                    %poles from expansion and mask
  allresid=pole2residue(allpole);
  [uhppole,uhpresid]=uhp(allpole,allresid);

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

end

[UHPpole,UHPresid]=uhp(ALLpole,ALLresid);
