%% bgden.m calculates the background H+ and He+ equatorial densities and
%% tube contents. This is done by assuming the densities are equal to 
%% values expected from the currently used polar model.
%%

function [hyeq,heeq,hytot,hetot]=bgden(lsh,kp)

% lsh = L-shell of field line of interest
% kp = 3-hourly Kp (scalar or same size as lsh)
% hyeq = equatorial H+ density (m^-3)
% heeq = equatorial He+ density (m^-3)
% hytot = H+ content of flux tube
% hetot = He+ content of flux tube

svm=ones(size(lsh));
re=6371;
dismat=re*lsh;
if prod(size(lsh))~=prod(size(kp))
 kp=kp*ones(size(lsh));
end
[hyeq,heeq]=polarmod(0*svm,dismat,kp);
vol=fltpr(lsh);
hytot=vol.*hyeq;
hetot=vol.*heeq;
