%% plsd42.m is designed to be called by plsd4.m. It calculates the
%% hydrogen ion density along the field lines from the base to the
%% maximum distance along the field line (or the field line apex if 
%% less). The profile is hybrid in the sense that it uses chemical
%% equilibrium from the base up to the point where the chemical 
%% loss/diffusion loss ratio is equal to the value inputted by 
%% the user, and from that point diffusive equilibrium is assumed. 
%% NOTE: 1) For values above the maximum, zero is returned. 
%%       2) Chemical equilibrium is assumed at r0. 

function hyden=plsd42(oxden,oxnt,hynt,dflm,g0,tmpi,tmpe,tmpt,tmpn,dflms,dfmx)

% oxden = oxygen ion density profiles along field lines (O+/m^3)
% oxnt = neutral oxygen density at points of interest (O/m^3)
% hynt = neutral hydrogen density at points of interest (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)
% dflms = distance along field line at which to swap from chemical to difussive control (km)
% dfmx = maximum distance along the field line to calculate H+ (scalar)(km)
%
% hyden = calculated hydrogen density values (H+/m^3)

%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
dflt=dflm(:,dvl);                       %length of the field line from the equator to r0
hyden=zeros(p1,dvl);                    %hydrogen density results matrix

% Determine initial chemical equilibrium H+ density at base
hyden(:,1)=oxden(:,1).*((2.5*hynt(:,1))./(2.2*oxnt(:,1))).*sqrt(tmpn(:,1)./tmpi(:,1));        %chemical equilibrium term at r0

z=1;                                    %parameter used control while loop
i=0;                                    %i'th section
dfl=zeros(p1,1);                        %distance of top of i'th section from base

while z
 
 i=i+1;                                 %increment i by 1
 ip=i+1;                                %i plus one
 
 % Determine the location of the top of the i'th section. This is done using
 % the distance along the field line from the base of field line.
 poi=dfl<dflt & dfl<dfmx;               %points of interest still above current dfl
 ddis=dflm(:,ip)-dflm(:,i);             %step size increment for i'th section 
 dfl(poi)=dflm(poi,ip);                 %new dfl values 
 
 % Determine if chemical equilibrium or diffusion to be used in i'th section
 poi1=poi & dflm(:,ip)<dflms;           %points where chemical equilibrium to be used
 poi2=poi & dflm(:,i)>=dflms;           %points where diffusive equilibrium to be used
 poi3=poi & ~poi1 & ~poi2;              %points where combined chemical and diffusive equilibrium to be used

 % Chemical equilibrium calculation
 if any(poi1)
  hyden(poi1,ip)=oxden(poi1,ip).*((2.5*hynt(poi1,ip))./(2.2*oxnt(poi1,ip))).*sqrt(tmpn(poi1,ip)./tmpi(poi1,ip));
 end 
 
 % Diffusion calculation                                                                                    
 if any(poi2)
  mb=(16*oxden(poi2,i)+hyden(poi2,i))./(oxden(poi2,i)+hyden(poi2,i));                     %mean ion density
  shhy=1e-3*(k*tmpi(poi2,i)./(amu*g0(poi2,i)))./(1-mb.*tmpe(poi2,i)./tmpt(poi2,i));       %hydrogen scale height (km)  
  hyden(poi2,ip)=hyden(poi2,i).*(tmpt(poi2,i)./tmpt(poi2,ip)).*exp(-ddis(poi2)./shhy);    %hydrogen density at top of section
 end

 % Combined calculation
 if any(poi3)
  hydence=oxden(poi3,ip).*((2.5*hynt(poi3,ip))./(2.2*oxnt(poi3,ip))).*sqrt(tmpn(poi3,ip)./tmpi(poi3,ip));     %chemical equilibrium calculation
  mb=(16*oxden(poi3,i)+hyden(poi3,i))./(oxden(poi3,i)+hyden(poi3,i));                                         %mean ion density
  shhy=1e-3*(k*tmpi(poi3,i)./(amu*g0(poi3,i)))./(1-mb.*tmpe(poi3,i)./tmpt(poi3,i));                           %hydrogen scale height (km)  
  hydende=hyden(poi3,i).*(tmpt(poi3,i)./tmpt(poi3,ip)).*exp(-ddis(poi3)./shhy);                               %hydrogen density at top of section
  wt=(dflms(poi3)-dflm(poi3,i))./ddis(poi3);                                              %weighting term
  hyden(poi3,ip)=wt.*hydence+(1-wt).*hydende;                                             %combined density profile
 end 
  
 % Check to see if any points of interest remain to be calculated. If not, break from loop
 if sum(dfl<dflt & dfl<dfmx)==0 | i==(dvl-1);
  z=0;
 end  
 
end 
