%% plsd43.m is designed to be called by plsd4.m. It calculates the
%% H+ flux at the distance diso along the field line. It uses the 
%% full diffusion eqn for the H+ flow from St-Maurice and Schunk, 
%% Planet. Space Sci. 1977. It is assumed that the O+ diffusion 
%% velocity is zero, and that the different ions have the same 
%% temperature. This is instead of simiplifed version given by 
%% Richards and Torr, 1985 JGR which plsd4 used previously.

function flux=plsd43(hyden,oxden,tmpe,tmpi,g0,dflm,diso)

% hyden = hydrogen ion density at points of interest (diffusive equilibrium) (H+/m^3)
% oxden = oxygen ion density at points of interest (O+/m^3)
% tmpe = electron temperature at points of interest (K)
% tmpi = ion temperature at points of interest (K)
% g0 = gravity at points of interest (m/s^2)
% dflm = distance along field lines of points of interest from r0 (km).
% diso = distance along the field line from r0 where flux calculated (km)
%
% flux = H+ flux when the diffusion loss/chemical loss ratio is one (H+/m^2) 

k=1.38e-23;                             %Boltzman constant (J/K) 
amu=1.66e-27;                           %atomic mass unit (kg)
re=6371;                                %Earth radi (km)
r0=re+350;                              %Base radial distance (km)

%Calculate the required parameters at diso
nop=size(dflm,1);                       %number of field lines
dfi=size(dflm,2);                       %number of points of interest per field line
poi=dflm<=(diso*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=diso-dflm(po1);                    %difference vector between diso value and lower point of interest

eden=oxden+hyden;                                 %e- density
m=(log(oxden(po2))-log(oxden(po1)))./dfv;
oxdeno=exp(m.*dfzf).*oxden(po1);                  %O+ density  
m=(log(hyden(po2))-log(hyden(po1)))./dfv;
hydeno=exp(m.*dfzf).*hyden(po1);                  %H+ density
dhyds=(m/1e3).*hydeno;                            %gradient of H+ density (H+/m^4)
m=(log(eden(po2))-log(eden(po1)))./dfv;
edeno=exp(m.*dfzf).*eden(po1);                    %e- density
deds=(m/1e3).*edeno;                              %gradient of e- density (e-/m^4)
m=(tmpi(po2)-tmpi(po1))./dfv;
tmpio=m.*dfzf+tmpi(po1);                          %ion temperature
dtids=m/1e3;                                      %gradient of ion temp (K/m)
m=(tmpe(po2)-tmpe(po1))./dfv;
tmpeo=m.*dfzf+tmpe(po1);                          %e- temperature
dteds=m/1e3;                                      %gradient of e- temp (K/m)
m=(g0(po2)-g0(po1))./dfv;
g0o=m.*dfzf+g0(po1);                              %%gravity (m/s^2)

%Determine estimates of alpha-alpha* and Delta using Fermi-Dirac funtion
%fit to values given by St-Maurice and Schunk for the O+/H+ case when
%the ion temperatures are equal.
nrat=log10(oxdeno./hydeno);                       %log10 of O+ to H+ ratio
acf=[1.2243 -2.3836 -0.3302 -0.4371];
dcf=[0.6370 -0.0001 -0.3301 -0.4344];
alpf=fedi(acf,nrat);
delf=fedi(dcf,nrat);

dfc=1e6*(k*tmpio.^(2.5))./(1.23*amu*oxdeno.*(1-delf));      %diffusion co-efficent (H+ m^2/s)

%Determine the values of the six terms in equation (19)
p1=dhyds./hydeno;
p2=amu*g0o./(k*tmpio);
p3=dtids./tmpio;
p4=(deds./edeno).*(tmpeo./tmpio);
p5=dteds./tmpio;
p6=-(oxdeno./edeno).*(dtids./tmpio).*alpf;

vel=-dfc.*(p1+p2+p3+p4+p5+p6);                         %velocity of H+ (m/s)
flux=vel.*hydeno;                                      %H+ flux (H+/m^2)

