%% plasipt16.m determines the He+ density down the given field lines. It can be
%% used in either the polar or mid/low latitude region. 
%% For regions below L=13 a combination of the diffusive equilibrium approach and 
%% the polar He+ profiles are used. For diffusive profiles with total content less
%% than polar profile content assume a polar density profile, for diffusive profiles
%% with total cotent great than five times the polar profile content assume a diffusive
%% profile and as a weight profile for diffusive contents in between. 
%% Above L=13, the He+ is always the polar profile. 
%% For points below L<13, plasipt16.m requires the H+ equatorial density, 
%% given by hyeq. For L>13 the hyeq value is not used.
%%
%% plasipt14.m does not calculate the conjugate half of the field line.
%% 
%% 23/11/2001 - modify all the profiles so that the final profiles equal to a weighted
%%              combination of the polar and diffusive profile, with the weighting
%%              based on total tube content.
%%  

function [heden,shhem]=plasipt16(oxden,hyden,heeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp)

% oxden = oxygen ion density profiles along field lines (O+/m^3)
% hyden = hydrogen ion density profiles along field lines (H+/m^3)
% heeq = equatorial helium ion density (H+/m^3) (value ignored for L>13) 
% lsh = L-shell of given field line (column vector)
% latmatm = magnetic latitude matrix of field lines (rads)
% dismat = radial displacement matrix of field lines (km)
% dflm = distance along field line of points on field line from base (km)
% g0 = gravity of points on field line (m/s^2)
% tmpi = ion temperature of points on field line (K)
% tmpe = electron temperature of points on field line (K)
% tmpt = total temperature of points on field line (K)
% kp = 3-hourly Kp (column vector)
%
% heden = calculated helium density values (He+/m^3)
% shhem = calculated helium scale height (km)

%Required constants
re=6371;                                %Earth's radius
k=1.38e-23;                             %Boltzman constant (J/K) 
amu=1.66e-27;                           %atomic mass unit (kg)

[p1,dvl]=size(oxden);                   %required size parameters

heden=zeros(p1,dvl);                    %H+ density results matrix
shhem=heden;                            %H+ scale height matrix
sv=ones(1,dvl);                         

poip=lsh>13;                            %polar He+ profiles
poi=~poip;                              %non-polar points

% Determine the polar profiles at all points of interest
[hypol,hepol]=polarmod(latmatm,dismat,kp);

% Non-polar diffusive equilibrium case
if any(poi)
 
 heden(poi,dvl)=heeq(poi);              %record non-polar equatorial densities
 
 for i=dvl:-1:2
  im=i-1;                               %i minus one
  if i==dvl 
   heden(poi,im)=heden(poi,i).*(tmpt(poi,i)./tmpt(poi,im));                                         %helium density at bottom of section
  else
   ddis=dflm(poi,i)-dflm(poi,im);                                                                   %length of i'th section
   mb=(16*oxden(poi,i)+4*heden(poi,i)+hyden(poi,i))./(oxden(poi,i)+heden(poi,i)+hyden(poi,i));      %mean ion density
   shhe=1e-3*(k*tmpi(poi,i)./(4*amu*g0(poi,i)))./(1-mb.*tmpe(poi,i)./(4*tmpt(poi,i)));              %helium scale height (km)  
   heden(poi,im)=heden(poi,i).*(tmpt(poi,i)./tmpt(poi,im)).*exp(ddis./shhe);                        %helium density at bottom of section
   shhem(poi,i)=shhe;                                                                               %record helium ion scale height (km)   
  end
 end 
end

%Determine polar diffusive equilibrium scale heights
if any(poip) 
 heden(poip,:)=hepol(poip,:);                                    %He+ polar densities
 for i=dvl:-1:2
  im=i-1;                                                        %i minus one
  if i~=dvl 
   mb=(16*oxden(poip,i)+4*heden(poip,i)+heden(poip,i))./(oxden(poip,i)+heden(poip,i)+heden(poip,i));     %mean ion density
   shhem(poip,i)=1e-3*(k*tmpi(poip,i)./(4*amu*g0(poip,i)))./(1-mb.*tmpe(poip,i)./(4*tmpt(poip,i)));      %helium scale height (km)  
  end
 end 
end

%%To determine the He+ density in the non-polar region combine the diffusive and
%%polar profiles using a weighting based on the polar and diffusive total content
%%of the tubes. They are joined using eqn(1) from Gallagher et al., JGR (2000).
if any(poi)
 hetotp=ftcont(hepol(poi,:),dismat(poi,:),dflm(poi,:),latmatm(poi,:));     %total polar He+ tube content
 hetot=ftcont(heden(poi,:),dismat(poi,:),dflm(poi,:),latmatm(poi,:));      %total diffusive He+ tube content
 rat=hetot./hetotp;                                                        %ratio of tube contents
 wt=0.5+0.5*tanh(3.4534*((rat-3)/2));                                      %weighting term
 heden(poi,:)=(1-wt*sv).*hepol(poi,:)+(wt*sv).*heden(poi,:);               %corrected He+ densities
end 

%Recalculate the He+ scale heights in non-polar region
if any(poi)
 for i=dvl:-1:2
  im=i-1;                                              %i minus one
  if i~=dvl 
   mb=(16*oxden(poi,i)+4*heden(poi,i)+heden(poi,i))./(oxden(poi,i)+heden(poi,i)+heden(poi,i));           %mean ion density
   shhem(poi,i)=1e-3*(k*tmpi(poi,i)./(4*amu*g0(poi,i)))./(1-mb.*tmpe(poi,i)./(4*tmpt(poi,i)));           %helium scale height (km)  
  end
 end 
end

%Calculate the helium scale height at r0
mb=(16*oxden(:,1)+4*heden(:,1)+heden(:,1))./(oxden(:,1)+heden(:,1)+heden(:,1));           %mean ion density 
shhem(:,1)=1e-3*(k*tmpi(:,1)./(4*amu*g0(:,1)))./(1-mb.*tmpe(:,1)./(4*tmpt(:,1)));         %helium scale height (km)
