%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             
% Name: LP_Vsc.m                                                              
% Author: Anders I Eriksson
% Modified: Kristoffer Hultgren
%
% Description: 
% Calculate spacecraft potential from plasma parameters                                                             
%                                                                         
%                                                                             
% Input:                                                                      
%   P.ne      -- [cm-3]
%   P.Te      -- [eV]
%   P.v       -- ram speed [km/s]
%   P.jf0     -- [A/m2]
%   vs0       -- seed value for Vsat [V]                                                                      
%                                                                             
% Output: vst -- [V] Satellite potential                                                                                                                                            
%                                                                             
%                                                                             
%                                                                             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vst = LP_Vsc(P,vs0)

global CO; % Global constants

vst = vs0;


v = (P.ni(1).*P.v(1) + P.ni(2).*P.v(2))./(P.ni(1) + P.ni(2));


Ie0 = P.ne*CO.qe*sqrt(P.Te/(2*pi*CO.me));
Ii0 = P.ne*CO.qe*v;
If0 = P.jf0;

if If0 + Ii0 < Ie0  % V < 0
    for iii=1:25 % Usually converges well within 5 steps
      vst = log(abs(1e3*v*(1-2*CO.qe*vst/(P.meff*CO.mp*1e6*v^2)) + P.jf0/(1e6*P.ne*CO.qe)));
      vst = vst - 0.5 * log(4*CO.qe*P.Te/(pi*CO.me)); %Plane geometry
      vst = P.Te * vst;
    end
end

if If0 + Ii0 > Ie0  % V > 0
    for iii=1:25 % Usually converges well within 5 steps
      vst = log(abs((1e6.*Ie0)./If0.*(1+CO.qe.*vst./P.Te) - 0.5.*(Ii0./If0.*(1-2.*CO.qe.*vst./(P.meff.*CO.mp.*1e6.*v^2)) + abs(Ii0./If0.*(1-2.*CO.qe.*vst/(P.meff.*CO.mp.*1e6.*v^2))))));
      vst = -P.Te .* vst;
    end
end