%% plsd62.m is designed to be called by plsd6.m. It calculates the
%% helium 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. 
%%
%% Unlike H+, He+ is primarily produced by photoionisation of neutral He
%% by sunlight. Since only considering altitudes above 350km assume no
%% change in solar intensity with altitude, therefore ionisation frequency
%% is integral value from ionfreq.m multipled by the neutral He density.
%% If the point of interest is in the solar shadow then no production 
%% occurs and the chemical equilibrium profile should be zero. However, for
%% GPID which is only concerned with altitudes above 350km it appears
%% that in the generation of He+ (at least) that using the chemical 
%% equilibrium profile generated with solar ionisation "shinning" during the
%% night still reproduces accurately enough the FLIP He+ profiles. Unlike 
%% with H, a smooth hybrid does not result when converting from chemical 
%% equilbrium straight to diffusive equilibrium, so instead use a weighting
%% from the swap height to swap height plus 200km.
%% 
%%
%% NOTE: 1) For values above the maximum, zero is returned. 
%%       2) Chemical equilibrium is assumed at r0.
%%

function heden=plsd62(oxden,n2nt,hent,dflm,g0,tmpi,tmpe,tmpt,tmpn,dflms,dise,f107,f107a,dfmx,tdiur)

% oxden = oxygen ion density profiles along field lines (O+/m^3)
% n2nt = neutral N2 density at points of interest (N2/m^3)
% hent = neutral helium 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)
% dise = distance along field line were chemical and diffusive loss are equal (km)
% f107 = current F10.7 solar flux (scalar or column vector)
% f107a = 81 day averaged F10.7 solar flux (scalar or column vector)
% dfmx = maximum distance along the field line to calculate He+ (scalar)(km)
% tdiur = diurinal parameter
%
% heden = calculated helium density values (He+/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
heden=zeros(p1,dvl);                    %helium density results matrix

%Determine the He+ ionisation frequency due to unattenuated solar flux.
if prod(size(f107))==1
 f107=f107*ones(p1,1);
 f107a=f107a*ones(p1,1);
end
qtot=zeros(p1,1);
for i=1:p1
 qtot(i)=ionfreq(f107(i),f107a(i),0);
end

% Determine initial chemical equilibrium H+ density at base
heden(:,1)=qtot.*hent(:,1)./(1.1e-15*n2nt(:,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, or a weighted combination of both.
 dflms2=(1-tdiur(:,i)).*(0.5*(dflms+dise))+tdiur(:,i).*dise;
 poi1=poi & dflm(:,ip)<=dflms;                         %points where chemical equilibrium to be used
 poi2=poi & dflm(:,ip)>=dflms2;                        %points where diffusive equilibrium to be used
 poi3=poi & dflm(:,ip)>dflms & dflm(:,ip)<dflms2;      %points where weighted combination to be used
 
 % Chemical equilibrium calculation
 if any(poi1)
  heden(poi1,ip)=qtot(poi1).*hent(poi1,ip)./(1.1e-15*n2nt(poi1,ip));
 end 
 
 % Diffusion calculation                                                                                    
 if any(poi2)
  mb=(16*oxden(poi2,i)+4*heden(poi2,i))./(oxden(poi2,i)+heden(poi2,i));                        %mean ion density
  shhe=1e-3*(k*tmpi(poi2,i)./(4*amu*g0(poi2,i)))./(1-mb.*tmpe(poi2,i)./(4*tmpt(poi2,i)));      %helium scale height (km)  
  heden(poi2,ip)=heden(poi2,i).*(tmpt(poi2,i)./tmpt(poi2,ip)).*exp(-ddis(poi2)./shhe);         %helium density at top of section
 end
 
 % Combined calculation
 if any(poi3)
  hece=qtot(poi3).*hent(poi3,ip)./(1.1e-15*n2nt(poi3,ip));                                     %chemical equilibrium densities
  mb=(16*oxden(poi3,i)+4*heden(poi3,i))./(oxden(poi3,i)+heden(poi3,i));                        %mean ion density
  shhe=1e-3*(k*tmpi(poi3,i)./(4*amu*g0(poi3,i)))./(1-mb.*tmpe(poi3,i)./(4*tmpt(poi3,i)));      %helium scale height (km)  
  hede=heden(poi3,i).*(tmpt(poi3,i)./tmpt(poi3,ip)).*exp(-ddis(poi3)./shhe);                   %helium density at top of section
  wt=(dflm(poi3,ip)-dflms(poi3))./(dflms2(poi3)-dflms(poi3));                                  %weighting term
  heden(poi3,ip)=(1-wt).*hece+wt.*hede;                                                        %weighted combined 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 

