% aa16case2bluewe input file for disprelelectronnew.
% This file is used to fit a simple pole distribution to the H2O+
% distribution function calculated from ICA data
%
% HG 2016-06-30

% ion mass is assumed to be in the variable mass and electron mass in me.
me     = 9.10938356e-31;
mH     = 1.672621898e-27;
mH2O   = 18*mH;
malpha = 4*mH;

mass_vec =[mH2O mH2O mH2O me]; % particle masses in kg
mass=mH2O;
mrat = mass/me;    % mass ratio
sqmr = sqrt(mrat); % square root of mass ratio


ALLpole =[];   
ALLresid=[];   

n_vec    =[4 2 2 5];                        % # of terms for Maxwell expansion
nb_vec   =[0 0 0 0];                        % # of terms for Butterworth      
op2      =[0.99986 0.00012 0.00002 mrat];   % densities 
vtbut_vec=[0 0 0 0]/sqmr;                   % Cut velocity in vte units
vdbut_vec=[0 0 0 0]/sqmr;                   % Butterworth center velocity
vt       =[0.34 0.3 0.5 sqmr]/sqmr;         % thermal velocity, vte=1
vd       =[0.0 3.85 8.1 0]/sqmr;            % drift velocity, vte=1



op2=(op2/sum(op2));
%%%% 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 i=1:n				% Poles for Maxwellian
    pol(2*(n-i)+1)=1.0/(2^i*prod(1:i));
  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 j=1:length(allpole)                  % return to variables with dimension
    ALLpole =[ALLpole allpole(j)*vt(jj)+vd(jj)];
    ALLresid=[ALLresid op2(jj)*allresid(j)];
  end
  
end

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