%% plsd4.m is designed to be called by the plasdyn#.m models. It estimates 
%% the H+ density profile at the bottom of the field lines by mixing the H+
%% diffusive profile projected down the field line with the hybrid diffusive/
%% chemical profile projected up from the base. The later 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. The two
%% profile mixing is controlled by the ratio of chemical loss/diffusion loss 
%% from the 1985 JGR papaer by Richards and Torr 1985. 
%% plsd4.m determines the the H+ flux at the point a distance diso along the
%% field line from the base. This is done using the full H+ drift velocity
%% from St-Maurice and Schunk, Planet Space Sci 1977 (eqn 19), assuming the
%% O+ drift velocity is zero.
%% NOTE - have decided to set diso equal to dise, but will retain the separate
%%        variables in case I choose to change this in the future.
%% NOTE - He+ is ignored when generating the upwards diffusive equilibrium
%%        component of the hybrid H+ profile. This is because O+ will be the
%%        dominant ion and He+ needs to be corrected in a similar manner to H+.
%%

function [hydenf,flux,diso,dise,hydenb,rat]=plsd4(hyden,oxden,shhy,shox,dismat,dflm,tmpi,tmpe,tmpt,tmpn,oxnt,hynt,g0,r0,rats)

% hyden = hydrogen ion density at points of interest (diffusive equilibrium) (H+/m^3)
% oxden = oxygen ion density at points of interest (O+/m^3)
% shhy = hydrogen 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)
% oxnt = neutral oxygen density at points of interest (O/m^3)
% hynt = neutral hydrogen density at points of interest (H/m^3)
% g0 = gravity along field line at points of interest (m/s^2)
% r0 = base radial distance of field lines (km) (scalar)
% rats = ratio at which to convert projected up profile from chemical to diffusive control (scalar)
%
% hydenf = modified hydrogen ion density profile (H+/m^3)
% flux = H+ flux when the diffusion loss/chemical loss ratio is one (H+/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)
% hydenb = H+ profile generated from projecting up from r0 (H+/m^3)
% rat = diffusion loss/chemical loss ratio

%% 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)=plsd41(dismat,dflm,tmpi,tmpn,oxden,shox,oxnt,shhy,r0,ratv(i)*sv);
end

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 H+ profile starting from r0
dfmx=5000;                                             %maximum height to calculate hybrid H+ 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
hydenb=plsd42(oxden,oxnt,hynt,dflm,g0,tmpi,tmpe,tmpt,tmpn,dflms,dfmx);

%% Combine the two H+ 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^(0.4)*rat^(1/5)
%% 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.4)*rat0(rat0>1e3).^(1/5);          %take rat^(1/5) for rat above 1000

hydenf=hyden;                                          %results matrix
wt=ones(size(tmpn));                                   %weighting matrix
wt(rat<0.01)=0;                                        %set chemical equilibrium values weighting to 0
hydenf(wt==0)=hydenb(wt==0);                           %record chemical equilibrium values

%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/5) 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=0.7;                                               
poic=log10(rat0)>(3-bet) & log10(rat0)<(3+bet);
if any(sum(poic))
 rt3=(-1/(30*bet))*log10(rat0(poic)).^2+(4/15+1/(5*bet))*log10(rat0(poic))+(1/5-bet/30-3/(10*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 
 hydenf(poi)=wt(poi).*hyden(poi)+(1-wt(poi)).*hydenb(poi);  %combined final H+ density
end

%% I have decided to add in one further correction. Calculate the H+/O+ ratio
%% at rat=1 and for all points at higher altitude. For points where the H+/O+
%% ratio has fallen below this value, reset H+ to maintain it. This to a first
%% approximation will enforce the fall (ie scale height) in H+ density is not
%% greater than that in O+, which is the greatest it can fall (p 666 Solar
%% Terrestrial Physics 555.514 S684).
nop=size(dflm,1);                                 %number of field lines
dfi=size(dflm,2);                                 %number of points of interest per field line
poi=dflm<=(dise*ones(1,dfi));                     %finds points in dflm matrix below or equal to z value
pois1=sum(poi,2);                                 %number of points in each row below or equal to z value 
pois2=pois1+1;                                    %pois1 plus 1
pois2(pois2==(dfi+1))=dfi;                        %reset required for points at 950km
po1=(pois1-1)*nop+(1:nop)';                       %index of pois1 points
po2=(pois2-1)*nop+(1:nop)';                       %index of pois2 points
dfv=dflm(po2)-dflm(po1);                          %difference vector between points of interest
dfzf=dise-dflm(po1);                              %difference vector between dise value and lower point of interest

m=(log(oxden(po2))-log(oxden(po1)))./dfv;
oxdeno=exp(m.*dfzf).*oxden(po1);                  %O+ density  

m=(log(hydenf(po2))-log(hydenf(po1)))./dfv;
hydenoh=exp(m.*dfzf).*hydenf(po1);  
horatb=(hydenoh./oxdeno)*ones(1,sz(2));                     %H+/O+ ratio at rat=0
hoc=dflm>(dise*ones(1,sz(2))) & (hydenf./oxden)<horatb;     %points to be corrected
hydenf(hoc)=oxden(hoc).*horatb(hoc);                        %corrected H+ values

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


