%% plasipt14.m determines the H+ 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 H+ 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 H+ is always the polar profile. 
%% For points below L<13, plasipt14.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 [hyden,shhym]=plasipt14(oxden,heden,hyeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp)

% oxden = oxygen ion density profiles along field lines (O+/m^3)
% heden = helium ion density profiles along field lines (He+/m^3)
% hyeq = equatorial hydrogen 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)
%
% hyden = calculated hydrogen density values (H+/m^3)
% shhym = calculated hydrogen 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

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

poip=lsh>13;                            %polar H+ 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)
 
 hyden(poi,dvl)=hyeq(poi);              %record non-polar equatorial densities
 
 for i=dvl:-1:2
  im=i-1;                               %i minus one
  if i==dvl 
   hyden(poi,im)=hyden(poi,i).*(tmpt(poi,i)./tmpt(poi,im));                                         %hydrogen 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
   shhy=1e-3*(k*tmpi(poi,i)./(amu*g0(poi,i)))./(1-mb.*tmpe(poi,i)./tmpt(poi,i));                    %hydrogen scale height (km)  
   hyden(poi,im)=hyden(poi,i).*(tmpt(poi,i)./tmpt(poi,im)).*exp(ddis./shhy);                        %hydrogen density at bottom of section
   shhym(poi,i)=shhy;                                                                               %record hydrogen ion scale height (km)   
  end
 end 
end

%Determine polar diffusive equilibrium scale heights
if any(poip) 
 hyden(poip,:)=hypol(poip,:);                               %H+ 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)+hyden(poip,i))./(oxden(poip,i)+heden(poip,i)+hyden(poip,i));     %mean ion density
   shhym(poip,i)=1e-3*(k*tmpi(poip,i)./(amu*g0(poip,i)))./(1-mb.*tmpe(poip,i)./tmpt(poip,i));            %hydrogen scale height (km)  
  end
 end 
end

%%To determine the H+ 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)
 hytotp=ftcont(hypol(poi,:),dismat(poi,:),dflm(poi,:),latmatm(poi,:));     %total polar H+ tube content
 hytot=ftcont(hyden(poi,:),dismat(poi,:),dflm(poi,:),latmatm(poi,:));      %total diffusive H+ tube content
 rat=hytot./hytotp;                                                        %ratio of tube contents
 wt=0.5+0.5*tanh(3.4534*((rat-3)/2));                                     %weighting term
 hyden(poi,:)=(1-wt*sv).*hypol(poi,:)+(wt*sv).*hyden(poi,:);               %corrected H+ densities
end 

%Recalculate the H+ scale heights in the 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)+hyden(poi,i))./(oxden(poi,i)+heden(poi,i)+hyden(poi,i));      %mean ion density
   shhym(poi,i)=1e-3*(k*tmpi(poi,i)./(amu*g0(poi,i)))./(1-mb.*tmpe(poi,i)./tmpt(poi,i));            %hydrogen scale height (km)  
  end
 end 
end

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

