%% plsd5.m calculates the required H+ production and loss arrays, from
%% which the change in H+ tube content can be determined. 
%% In calculating the downwards diffusion flux, any positive terms are 
%% reset to zero, since positive values are assumed to correspond to the 
%% nett upwards H+ flow. The largest allowed downward diffusive flux is 
%% set to -3.5e12 H+/m^2 s^-1, which is larger than the value noted by 
%% Park 1970 based on whistler measurements (note that this was for L=4 
%% at 1000km) of 1.5e12, and is roughly the upper upwards flux observed
%% by Park and others.
%% Note that this correction is done before the value is modifed due to
%% the cross sectional area of the given flux tube of interest.
%% A maximum value is placed on the nett flux, given by eqn 4 from 
%% Richards and Torr 1985. The maximum value derived from this equation 
%% is modifed depending on the F107 (this is based on comparisons to
%% FLIP results and looking at plots in Richards and Torr comparing 
%% eqn 4 to the full FLIP model). At low F107 the flux is reduced, while
%% at high F107 it is increased.
%% NOTE - the Richards and Torr equation was compared to the flux through
%%        3000km (total H+ produced in flux tube divided by cross sectional
%%        area at 3000km) then normalised to 200km. Thus the H+ fluxes are 
%%        effectively equal to the total production in the flux tube, since
%%        base cross-sectional areas all roughly 1m^2, and this doesn't
%%        change by much moving up to 350km used by GPID.
%%

function [plsf,hfluxt,pdfi,lsfi,plsfmx,hychem]=plsd5(lsh,latmatm,dflm,tmpi,tmpn,hyden,oxden,oxnt,hynt,shox,csa,hflux,diso,dise,kp,f107)

% lsh = L-shells of interest
% latmatm = magnetic latitude of profile points (rads)
% dflm = distance along field line from r0 of profile points (km)
% tmpi = ion temperature at points of interest (K)
% tmpn = neutral temperature at points of interest (K)
% hyden = hydrogen ion density at points of interest (H+/m^3)
% oxden = oxygen ion density at points of interest (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)
% shox = O+ scale height (km)
% csa = cross sectional area of each point on the field line (m^2)
% hflux = 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 chemical loss is equal to diffusive loss (km)
% kp = 3-hourly Kp solar activity (scalar)
% f107 = F10.7 solar flux (scalar)
%
% plsf = net upwards H+ flux (H+/s)
% hfluxt = total maximum downward H+ diffusive flux for flux tube (H+/s)
% pdfi = sectional H+ flux production along field line (H+/s)
% lsfi = sectional H+ flux loss along field line (H+/s)
% plsfmx = maximum upwards flux from Richards and Torr eqn 5 (H+/s)
% hychem = net upwards H+ flux before inclusion of Richards and Torr limit (H+/s)

re=6371;                                %Earth's radius
r0=re+350;                              %Base radial distance

%% Calculate the H+ chemical production and loss occuring along the field line
dvl=size(dflm,2);
dvl2=dvl-1;
ddis=1e3*(dflm(:,2:dvl)-dflm(:,1:dvl2));                         %width of each section (m)

lsfp=2.2e-17*sqrt(tmpi).*oxnt.*hyden;                                                %loss at each point (H+/m^3/s)
lsfi=0.5*(lsfp(:,(1:dvl2)).*csa(:,(1:dvl2))+lsfp(:,(2:dvl)).*csa(:,(2:dvl))).*ddis;  %loss term in each section (H+/s)
pdfp=2.5e-17*sqrt(tmpn).*hynt.*oxden;                                                %production at each point (H+/m^3/s)
pdfi=0.5*(pdfp(:,(1:dvl2)).*csa(:,(1:dvl2))+pdfp(:,(2:dvl)).*csa(:,(2:dvl))).*ddis;  %production term in each section (H+/s)
plsf=sum(pdfi-lsfi,2);                                                               %calculate net upward flux

%Record nett production and loss before using upper limit
hychem=plsf;

%%Determine the maximum allowed nett flux
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 (O+/m^3) 
m=(log(hynt(po2))-log(hynt(po1)))./dfv;
hynto=exp(m.*dfzf).*hynt(po1);                    %H density  (H+/m^3)
m=(shox(po2)-shox(po1))./dfv;
shoxo=1e3*(m.*dfzf+shox(po1));                    %O+ scale height (m)
m=(tmpn(po2)-tmpn(po1))./dfv;
tmpno=m.*dfzf+tmpn(po1);                          %neutral temperature (K)

plsfmx=2.5e-17*sqrt(tmpno).*oxdeno.*hynto.*shoxo;      %maximum H+ nett flux

%Modify the maximum flux using the F107. Use a linear weighting term.
%Further, scale the maximum flux by Kp as well, using the weighting 
%function given by lshsf
wt=0.7+(0.1/30)*(f107-70);
if kp<=1.5
 kpsf=1;
else
 kpsf=(0.33*(kp-1.5)*lshsf(lsh))+1;
end
plsfmx=wt.*plsfmx./kpsf;
plsf(plsf>plsfmx)=plsfmx(plsf>plsfmx);                 %reset required maximum flux

%%Determine the total diffusive flux loss along the field line. Note that each flux
%%is considered to be zero if negative ie ignore positive since we are considering
%%that any positive values corresponding to the upward flow of the nett H+ production.
%%The maximum allowed flux is -3.5e12. Note that this is the maximum flux per m^2, ie 
%%before effect of the cross-sectional area taken into account for each flux tube.
hflux(hflux>0)=0;                                      %reset flux values above zero to zero
hflux(hflux<-3.5e12)=-3.5e12;

latmf=dis2th(diso,latmatm(:,1),r0);                    %magnetic latitude of flux point
disf=lshpar(lsh,latmf,-ones(size(latmf)));             %radial distance to flux point
bo=30.4e-6;                                            %equatorial magnetic field (T)
lshs=3;                                                %standard L-shell
ths=acos(sqrt((1/lshs)*(r0/re)));                      %latitude of standard point (rads)
fls=bo*(re/r0)^3*sqrt(1+3*sin(ths)^2);                 %standard magnetic flux (T.m^2)
bf=bo*(re./disf).^3.*sqrt(1+3*sin(latmf).^2);          %magnetic field at flux point (T)
hfluxt=hflux.*fls./bf;                                 %total H+ flux at flux point (H+/s)
