%% plsd3.m is designed to be called by the plasdyn#.m models. It calculates
%% the helium ion density along the field lines starting at the equator and
%% project down the field line using diffusive equilibrium. 
%%

function [heden,hedenc,shhem,shhemc]=plsd3(lsh,oxden,oxdenc,hyden,hydenc,heeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc,tdiur,tdiurc,f107)

% lsh = L-shells of interest (column vector)
% 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)
% hyden = hydrogen ion density profiles along field lines (H+/m^3)
% hydenc = conjugate hydrogen ion density profiles along field lines (H+/m^3)
% heeq = helium equatorial ion density (He+/m^3)  (column vector)
% 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)
% tdiur = diurinal parameters
% tdiurc = conjugate diurinal parameters
% f107 = F10.7 solar flux  (column vector) 
%
% heden = calculated helium density values (He+/m^3)
% hedenc = calculated conjugate helium density values (He+/m^3)
% shhem = calculated helium scale height (km)
% shhemc = calculated conjugate 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);                    %helium density values
shhem=heden;                            %helium scale heights
shhemc=heden;                           %conjugate helium scale heights
heden(:,dvl)=heeq;                      %record equatorial helium density
hedenc=heden;                           %conjugate helium 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 
  heden(:,im)=heden(:,i).*(tmpt(:,i)./tmpt(:,im));                                   %helium density at bottom of section
 else 
  mb=(16*oxden(:,i)+4*heden(:,i)+hyden(:,i))./(oxden(:,i)+heden(:,i)+hyden(:,i));    %mean ion density
  shhe=1e-3*(k*tmpi(:,i)./(4*amu*g0(:,i)))./(1-mb.*tmpe(:,i)./(4*tmpt(:,i)));        %helium scale height (km)  
  shhe=shhe.*hesf(f107,tdiur(:,i),lsh);                                              %scale the He+ using the He+ scaling factor
  heden(:,im)=heden(:,i).*(tmpt(:,i)./tmpt(:,im)).*exp(ddis./shhe);                  %helium density at bottom of section
  shhem(:,i)=shhe;                                                                   %record helium ion scale height (km)
 end
 
 %%Second consider the conjugate point of interests
 if i==dvl 
  hedenc(:,im)=hedenc(:,i).*(tmptc(:,i)./tmptc(:,im));                                    %helium density at bottom of section
 else 
  mb=(16*oxdenc(:,i)+4*hedenc(:,i)+hydenc(:,i))./(oxdenc(:,i)+hedenc(:,i)+hydenc(:,i));   %mean ion density
  shhe=1e-3*(k*tmpic(:,i)./(4*amu*g0(:,i)))./(1-mb.*tmpec(:,i)./(4*tmptc(:,i)));          %helium scale height (km)  
  shhe=shhe.*hesf(f107,tdiurc(:,i),lsh);                                                  %scale the He+ using the He+ scaling factor
  hedenc(:,im)=hedenc(:,i).*(tmptc(:,i)./tmptc(:,im)).*exp(ddis./shhe);                   %helium density at bottom of section
  shhemc(:,i)=shhe;                                                                       %record helium ion scale height (km)
 end
  
end

%Calculate the helium scale height at r0
mb=(16*oxden(:,1)+4*heden(:,1)+hyden(:,1))./(oxden(:,1)+heden(:,1)+hyden(:,1));           %mean ion density 
shhe=1e-3*(k*tmpi(:,1)./(4*amu*g0(:,1)))./(1-mb.*tmpe(:,1)./(4*tmpt(:,1)));               %helium scale height (km)
shhem(:,1)=shhe.*hesf(f107,tdiur(:,1),lsh);                                               %scale the He+ using the He+ scaling factor
mb=(16*oxdenc(:,1)+4*hedenc(:,1)+hydenc(:,1))./(oxdenc(:,1)+hedenc(:,1)+hydenc(:,1));     %mean ion density 
shhe=1e-3*(k*tmpic(:,1)./(4*amu*g0(:,1)))./(1-mb.*tmpec(:,1)./(4*tmptc(:,1)));            %conjugate helium scale height (km)
shhemc(:,1)=shhe.*hesf(f107,tdiurc(:,1),lsh);                                             %scale the He+ using the He+ scaling factor

