%% plsd7.m calculates the required He+ production and loss arrays, from
%% which the change in He+ tube content can be determined. It is based on
%% plsd5.m, which claculates the corresponding H+ production.
%% In calculating the downwards diffusion flux, any positive terms are 
%% reset to zero, since positive values are assumed to correspond to the 
%% nett upwards He+ flow. 
%% Because He+ is produced by photoionisation, need to determine if given
%% point along a field line is in the Earth's shadow. For those that are
%% the ionisation frequency is set to zero.
%%

function [plsf,hefluxt,pdfi,lsfi]=plsd7(lsh,dismat,latmatm,dflm,csa,lonmat,latmat,timmat,heden,n2nt,hent,g0,tmpn,heflux,diso,dise,kp,f107,f107a,year,month,day)

% lsh = L-shells of interest
% dismat = radial distance of points on field line from base (km)
% latmatm = magnetic latitude of profile points (rads)
% dflm = distance along field line from r0 of profile points (km)
% csa = cross sectional area of each point on the field line (m^2)
% lonmat = geographic longitude of each point on the field line (rads)
% latmat = geographic latitude of each point on the field line (rads)
% timmat = univerisal time (fractional hours) (scalar)
% heden = helium ion density at points of interest (He+/m^3)
% n2nt = neutral N2 density at points of interest (n2/m^3)
% hent = neutral helium density at points of interest (He/m^3)
% g0 = gravity of points on field line (m/s^2)
% tmpn = neutral temperature of points on field line (K)
% heflux = 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 chemical loss is equal to diffusive loss (km)
% kp = 3-hourly Kp solar activity (scalar)
% f107 = F10.7 solar flux (scalar)
% f107a = 81 day average F10.7 solar flux (scalar)
% year = starting year (scalar)
% month = starting month (scalar)
% day = starting day (scalar)
%
% plsf = net upwards He+ flux (He+/s)
% hefluxt = total maximum downward He+ diffusive flux for flux tube (He+/s)
% pdfi = sectional He+ flux production along field line (He+/s)
% lsfi = sectional He+ flux loss along field line (He+/s)

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

%Determine the He+ ionisation frequency due to unattenuated solar flux.
qtot=ionfreq(f107,f107a,0);

%Calculate if each point is in the Earth's "shadow".
svm=ones(size(lonmat));
timl=timut2l(timmat*svm,lonmat);
pois=eshadow(day,month,year,timl,latmat,dismat);

%Calculate N2 scale height at dise for each field line
tmpnde=intfl(dflm,tmpn,dise,1);         %neutral temperature at dise (K)
g0de=intfl(dflm,g0,dise,1);             %gravity at dise (m.s^-2)
n2sh=1e-3*k*tmpnde./(28*amu*g0de);      %N2 sacle height at dise (km)

%Calculate which points are used below the dise location for production and loss calculation
p1=[80 2];
p2=[200 0];
m=(p2(2)-p1(2))/(p2(1)-p1(1));
ssf=p1(2)+m*(f107-p1(1));                              %Solar scaling factor due to F107
if kp<=1.5
 kpsf=2*ones(size(lsh));                               %Kp scaling factor 
else
 kpsf=0.5*(kp-1.5)*lshsf(lsh)+2;
end
poilos=dflm>((dise-(kpsf.*n2sh))*ones(1,size(dflm,2)));
poipro=dflm>((dise-(ssf.*n2sh))*ones(1,size(dflm,2)));

%% Calculate the He+ 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=1.1e-15*n2nt.*heden.*poilos;                                                         %loss at each point (He+/m^3/s)
lsfi=0.5*(lsfp(:,(1:dvl2)).*csa(:,(1:dvl2))+lsfp(:,(2:dvl)).*csa(:,(2:dvl))).*ddis;       %loss term in each section (He+/s)
pdfp=qtot*hent.*(~pois).*poipro;                                                          %production at each point (He+/m^3/s)
pdfi=0.5*(pdfp(:,(1:dvl2)).*csa(:,(1:dvl2))+pdfp(:,(2:dvl)).*csa(:,(2:dvl))).*ddis;       %production term in each section (He+/s)
plsf=sum(pdfi-lsfi,2);                                                                    %calculate net upward flux

%After looking at refilling comparisons (flipcp2.m) is appears that a correction factor at low L 
%is needed to increase He+ diffusive flux.
cv=[23 -95 98];
lscf=cv(1)*lsh.^2+cv(2)*lsh+cv(3);
lscf0=cv(1)*2.^2+cv(2)*2+cv(3);
lsc=1+lscf-lscf0;
lsc(lsh>2)=1;
heflux=lsc.*heflux;

%%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 He+ production.
heflux(heflux>0)=0;                                    %reset flux values above zero to zero
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)
hefluxt=heflux.*fls./bf;                               %total He+ flux at flux point (He+/s)

