%% plsd2.m is designed to be called by the plasdyn#.m models. It calculates
%% the hydrogen ion density along the field lines starting at the equator and
%% project down the field line using diffusive equilibrium. 
%%
%% plsd2.m is based on the earlier plsd9.m, but now includes He+.
%%

function [hyden,hydenc,shhym,shhymc]=plsd2(oxden,oxdenc,heden,hedenc,hyeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc)

% oxden = oxygen ion density profiles along field lines (initial column has value at base of field line) (O+/m^3)
% oxdenc = conjugate oxygen ion density profiles along field lines (initial column has value at base of field line) (O+/m^3)
% heden = helium ion density profiles along field lines (He+/m^3)
% hedenc = conjugate helium ion density profiles along field lines (He+/m^3)
% hyeq = hydrogen equatorial ion density (H+/m^3)
% 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)
% tmpn = neutral temperature of points on field line (K)
% tmpic = conjugate ion temperature of points on field line (K)
% tmpec = conjugate electron temperature of points on field line (K)
% tmptc = conjugate total temperature of points on field line (K)
% tmpnc = neutral temperature of points on field line (K)
%
% hyden = calculated hydrogen density values (H+/m^3)
% hydenc = calculated conjugate hydrogen density values (H+/m^3)
% shhym = calculated hydrogen scale height (km)
% shhymc = calculated conjugate 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);                    %hydrogen density values
shhym=hyden;                            %hydrogen scale heights
shhymc=hyden;                           %conjugate hydrogen scale heights
hyden(:,dvl)=hyeq;                      %record equatorial hydrogen density
hydenc=hyden;                           %conjugate hydrogen density values

for i=dvl:-1:2
 
 im=i-1;                                %i minus one
 ddis=dflm(:,i)-dflm(:,im);             %length of i'th section
  
 %%First consider the point of interests
 if i==dvl 
  hyden(:,im)=hyden(:,i).*(tmpt(:,i)./tmpt(:,im));                                   %hydrogen density at bottom of section
 else 
  mb=(16*oxden(:,i)+4*heden(:,i)+hyden(:,i))./(oxden(:,i)+heden(:,i)+hyden(:,i));    %mean ion density
  shhy=1e-3*(k*tmpi(:,i)./(amu*g0(:,i)))./(1-mb.*tmpe(:,i)./tmpt(:,i));              %hydrogen scale height (km)  
  hyden(:,im)=hyden(:,i).*(tmpt(:,i)./tmpt(:,im)).*exp(ddis./shhy);                  %hydrogen density at bottom of section
  shhym(:,i)=shhy;                                                                   %record hydrogen ion scale height (km)
 end
 
 %%Second consider the conjugate point of interests
 if i==dvl 
  hydenc(:,im)=hydenc(:,i).*(tmptc(:,i)./tmptc(:,im));                                    %hydrogen density at bottom of section
 else 
  mb=(16*oxdenc(:,i)+4*hedenc(:,i)+hydenc(:,i))./(oxdenc(:,i)+hedenc(:,i)+hydenc(:,i));   %mean ion density
  shhy=1e-3*(k*tmpic(:,i)./(amu*g0(:,i)))./(1-mb.*tmpec(:,i)./tmptc(:,i));                %hydrogen scale height (km)  
  hydenc(:,im)=hydenc(:,i).*(tmptc(:,i)./tmptc(:,im)).*exp(ddis./shhy);                   %hydrogen density at bottom of section
  shhymc(:,i)=shhy;                                                                       %record hydrogen ion scale height (km)
 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)
mb=(16*oxdenc(:,1)+4*hedenc(:,1)+hydenc(:,1))./(oxdenc(:,1)+hedenc(:,1)+hydenc(:,1));     %mean ion density 
shhymc(:,1)=1e-3*(k*tmpic(:,1)./(amu*g0(:,1)))./(1-mb.*tmpec(:,1)./tmptc(:,1));           %conjugate hydrogen scale height (km)


