%% plsd6.m is designed to be called by the plasdyn#.m models. Based on plsd4.m it 
%% calculates the He+ profiles in the topside ionosphere. Following the same rational
%% as plsd4.m, in plsd6.m the H+ is ignored when generating the upward He+ diffusive
%% equilibrium component of the hybrid profile.
%%

function [hedenf,flux,diso,dise,hedenb,rat]=plsd6(heden,oxden,shhe,shox,dismat,dflm,tmpi,tmpe,tmpt,tmpn,n2nt,hent,g0,tdiur,f107,f107a,r0,rats)

% heden = helium ion density at points of interest (diffusive equilibrium) (He+/m^3)
% oxden = oxygen ion density at points of interest (O+/m^3)
% shhe = helium ion scale heights at points of interest (km)
% shox = oxygen ion scale heights at points of interest (km)
% dismat = distance from center of Earth of points of interest (km).
% dflm = distance along field lines of points of interest from r0 (km).
% tmpi = ion temperature at points of interest (K)
% tmpe = electron temperature at points of interest (K)
% tmpt = ion plus electron temperature (K)
% tmpn = neutral temperature at points of interest (K)
% n2nt = neutral nitrogen molecular (N2) density at points of interest (N2/m^3)
% hent = neutral helium density at points of interest (He/m^3)
% g0 = gravity along field line at points of interest (m/s^2)
% tdiur = diurinal parameter
% f107 = current F10.7 solar flux (scalar or column vector)
% f107a = 81 day averaged F10.7 solar flux (scalar or column vector)
% r0 = base radial distance of field lines (km) (scalar)
% rats = ratio at which to convert projected up profile from chemical to diffusive control (scalar)
%
% hedenf = modified helium ion density profile (He+/m^3)
% flux = He+ flux when the diffusion loss/chemical loss ratio is one (He+/m^2) 
% diso = distance along the field line from r0 where flux calculated (km)
% dise = distance along the field line from r0 where diffusion loss/chemical loss are equal (km)
% hedenb = He+ profile generated from projecting up from r0 (He+/m^3)
% rat = diffusion loss/chemical loss ratio

%Multiple original rats by 0.1
rats=0.1*rats;

%% First need to determine the distance along the field lines at which the ratio of 
%% of diffusion loss/chemical loss is 10 and 0.1. Using linear interpolation,
%% determine the ratio at all other points of interest.
ratv=[0.1 10];                                         %required ratio values
sz=size(dismat);                                       %size of input
z=zeros(sz(1),2);                                      %results matrix
sv=ones(sz(1),1);

for i=1:2                                              %estimate required distances
 z(:,i)=plsd61(dismat,dflm,tmpi,tmpn,oxden,shox,n2nt,shhe,r0,ratv(i)*sv);
end

%Check for z(:,2)-z(:,1)<50 (can occur at times). Reset these values so maximum z(:,1) is z(:,2)-50
poiz=z(:,2)-z(:,1)<50;
z(poiz,1)=z(poiz,2)-50;

mr=(z(:,2)-z(:,1))./((log10(ratv(2))-log10(ratv(1)))*sv);   %ratio gradient term
mr2=mr;                                                     %copy of mr
mr2(mr2==0)=-1;                                             %reset terms with same z values to 'marker' value (required for very low L-shell terms)
sv2=ones(1,sz(2));
rat=ratv(1).*10.^((dflm-z(:,1)*sv2)./(mr2*sv2));       %diffusion loss/chemical loss ratios
rat(mr2==-1,:)=0.01;                                   %reset required terms (occurs only at very low L-shell ie L<1.1)

%Restrict extreme low rat values
rat(rat<0.001)=0.001;

%% Determine distance along field line where ratio is one
dise=z(:,1)-mr.*log10(ratv(1));                                

%% Calculate the hybrid He+ profile starting from r0
dfmx=1e6;                                              %maximum height to calculate hybrid He+ up to
dflms=z(:,1)+mr.*log10(rats./ratv(1));                 %distance along field line at which to swap lower profile from chemical to diffusive control

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

%% Combine the two He+ profiles using the rat ratio values. After comparison with 
%% some FLIP profiles, I have found that better agreement can be obtained for values
%% with rat>1 by using rat^(1/3). Further, for values above rat=1000 use 10^(1/3-1/10)*rat^(1/10)
%% where 10^(0.4) is so two profiles agree at rat=1000.
rat0=rat;                                                   %copy of orignal ratios
rat(rat0>1 & rat<=1e3)=rat0(rat0>1 & rat<=1e3).^(1/3);      %take rat^(1/3) for rat between 1 and 1000
rat(rat0>1e3)=10^(0.7)*rat0(rat0>1e3).^(1/10);              %take rat^(1/10) for rat above 1000

hedenf=heden;                                          %results matrix
wt=ones(size(tmpn));                                   %weighting matrix

%Correction factor to smoothly join r with r^(1/3) in region around rat=1. Assume that
%this is given by a parabola fitted between log10(rat0)=-alp and log10(rat0)=alp were
%rat0 are the original ratios.
alp=0.7;                                               
poic=log10(rat0)>-alp & log10(rat0)<alp;
if any(sum(poic))
 rt2=(-1/(6*alp))*log10(rat0(poic)).^2+(2/3)*log10(rat0(poic))-alp/6;
 rat(poic)=10.^rt2;
end 

%Correction factor to smoothly join r^(1/3) with r^(1/10) in region around rat=1000. Assume 
%that this is given by a parabola fitted between log10(rat0)=-bet and log10(rat0)=bet were
%rat0 are the original ratios.
bet=1.4;                                               
poic=log10(rat0)>(3-bet) & log10(rat0)<(3+bet);
if any(sum(poic))
 rt3=(-7/(120*bet))*log10(rat0(poic)).^2+(13/60+7/(20*bet))*log10(rat0(poic))+(7/20-bet*(7/120)-21/(40*bet));
 rat(poic)=10.^rt3;
end 

%Combine the two profiles
poi=rat>0.01 & rat<1e3;                                     %terms to be combined
if any(sum(poi))
 wt(poi)=rat(poi)./(rat(poi)+1);                            %weighting term for poi 
 hedenf(poi)=wt(poi).*heden(poi)+(1-wt(poi)).*hedenb(poi);  %combined final H+ density
end

%Values below rat=0.01 are assumed to be in chemical equilibrium
hedenf(rat<=0.01)=hedenb(rat<=0.01);                        %record chemical equilibrium values

%% Determine the He+ flux at diso in the final corrected He+ profile.
diso=dise;                                                  %set diso equal to dise
flux=plsd63(hedenf,oxden,tmpe,tmpi,g0,dflm,diso);
