%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                           
% Name: LP_Curri.m                          
% Author: Anders I. Eriksson
%
% Description:                                                                        
%
%  Model uses Kanal's (1962) results as cited by Whipple 
% (1965, eqn 3.9); these are identical to Hoegy & Brace 
% (1999) expressions. Medicus (1962) expressions should
% also be identical in the limit of large sheath radius,
% though this remains to check. 
%
% These should always apply for repulsive potential, but
% for attractive potential they are only correct for OML.
% By supplying the Medicus (1962) results with a relation
% for the sheath radius, it should be possible to get a
% better model for the case of a spherical sheath. The
% more realistic case of a wake requires expressions
% fitted to simulation results.
%
% Input:                                    
%                                           
%  vp: probe potential [V]
%   n: ion density [cm-3]
%   T: ion temperature [eV]
%   z: ion charge [elementary charge]
%   m: ion mass [proton masses]
%   v: ram speed [km/s]
%   dn: dn/n noise [%]
%                                  
% Output:                                 
%                                         
%   ii: current [A]                                   
%                                         
% Notes:                                  
%                                         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ii = LP_Curri(vp,n,T,z,m,v,dn)

global IN CO;

q = z*CO.qe;

n = n*1e6;
v = v*1e3;
m = m*CO.mp;

% Random current:
ii0 = IN.probe_radius^2*n*(1+dn/100)*q*sqrt(8*pi*CO.qe*T/m);

att = find(q*vp < 0);   % Attracted ions
rep = find(q*vp >= 0);  % Repelled ions

% Initialize return:
ii = zeros(size(vp));

vt = sqrt(2*CO.qe*T/m);  % Thermal speed [m/s]
r = v/vt;

% Repelled ions:
eta = -z*vp(rep)/T;
sqe = sqrt(-eta);
tmp1 = (eta + r^2 + 1/2) * (sqrt(pi)/2);
tmp2 = erf(sqe + r) - erf(sqe - r);
tmp3 = (sqe + r) .* exp( -(sqe - r).^2 )/2;
tmp4 = (sqe - r) .* exp( -(sqe + r).^2 )/2;
ii(rep) = (1/(2 * r)) * (tmp1 .* tmp2 + tmp3 - tmp4);

% Attracted ions:
eta = -z*vp(att)/T;
sqe = sqrt(eta);
tmp1 = exp(-r^2)/2;
tmp2 = sqrt(pi) * (eta + r^2 +1/2) * erf(r) / (2*r);
ii(att) = (tmp1 + tmp2);

ii = -ii0 .* ii;

