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

ALLpole =[];
ALLresid=[];

n_vec =[2 4 2 3]; % # of terms for Maxwell expansion
nb_vec =[0 0 0 0]; % # of terms for Butterworth
op2 =[0.801864 0.045007 0.150910 0.002219]; % densities
vtbut_vec=[0 0 0 0]; % Cut velocity in cs units
vdbut_vec=[0 0.0 0.239 0]; % Butterworth center velocity
vt =[0.0359 0.275 0.27 0.478]; % thermal velocity, cs=1
vd =[0.0 2.63 3.6 5.5]; % drift velocity, cs=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