%% plasdyn3.m is an update of plasdyn8.m that calculates the oxygen, hydrogen and helium
%% ion densities profiles at the selected heights along the given L-shells. It is designed
%% to be called by plasimod1.m and should probably not be called otherwise as it requires a
%% specific input format.
%%  
%% The L-shell is assumed to have its base at 350km altitude. It then calculates the
%% H+ and He+ production and loss along the field line to determine the net gain or loss in 
%% the respective ions tube content. All output is in the format such that each row corresponds
%% to a seperately inputted field line of interest, except the new total tube contents
%% and new equatorial densities are outputted with the same size as the input.
%% 
%% NOTE 1) the co-ordinates, temperatures and densities parameters require global 
%%         variables (ie they have to be defined and calculated by the routine
%%         that calls plasdyn8.m).
%%      

function [hytotn,hetotn,oxeqf,hyeqf,heeqf,hyden,heden,oxden,hydenc,hedenc,oxdenc]=plasdyn8(poi,hytot,hetot,t0d,t0n,g0d,g0n,t0dc,t0nc,g0dc,g0nc,sl,bh,dhr,lsh,timi,timut,dno,year,tint,lshcp,f107,f107a,kp,hes)

% poi = logic vector to select required field lines from global variables
% hytot = initial H+ total tube content of required field lines
% hetot = initial He+ total tube content of required field lines
% t0d = day base temperature (K)
% t0n = night base temperature (K)
% g0d = day gradient temperature (K/km)
% g0n = night gradient temperature (K/km)
% t0dc = conjugate day base temperature (K)
% t0nc = conjugate night base temperature (K)
% g0dc = conjugate day gradient temperature (K/km)
% g0nc = conjugate night gradient temperature (K/km)
% sl = sL parameter vector
% bh = bh parameter vector
% dhr = dhr parameter vector
% lsh = L-shell vector for one magnetic longitude
% timi = time index to be considered
% timut = UT (scalar)
% dno = day number (scalar)
% year = year (scalar)
% tint = time interval to calculate change in tube content (fractional hours) (scalar)
% lshcp = maximum L-shell for calculation of profiles (scalar)
% f107 = F107 solar flux (scalar)
% f107a = 81 day averaged F107 solar flux (scalar)
% kp = Kp (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
%
% The following GLOBAL VARIABLES are required data for field line co-ordinates
% LON350 = geographic longitude at 350km (rads)
% LAT350 = geographic latitude at 350km (rads)
% LON350C = conjugate geographic longitude at 350km (rads)
% LAT350C = conjugate geographic latitude at 350km (rads)
% LATMATM = magnetic latitude (rads)
% DISMAT = radial distance (km)
% DFLM = distance along field line from base
% LSHV = L-shell of field line (vector)
%
% The following GLOBAL VARIABLES (+ (sometimes) conjugates ) are required data for electron density
% EDEN350 = electron density at 350km (m^-3)
% HMF2 = HmF2 (km)
% NMF2 = NmF2 (m^-3)
% BETA, ETA, DELTA, ZETA = Bent model parameters
%
% The following GLOBAL VARIABLES (+ (sometimes) conjugates ) are required data for the neutral parameters
% OXNT350 = neutral oxygen density at 350km
% HYNT350 = neutral hydrogen density at 350km
% HENT350 = neutral helium density at 350km
% N2NT350 = neutral molecular nitrogen density at 350km
% TMPN350 = neutral temperature at 350km
% TMPNSP = neutral exospheric temperature
% TMPNSP = neutral temperature profile shaping parameter
% 
% hytotn = new total hydrogen ion content of flux tube 
% hetotn = new total helium ion content of flux tube 
% oxeqf = final equatorial oxygen ion density of flux tube (m^-3)
% hyeqf = final equatorial hydrogen ion density of flux tube (m^-3)
% heeqf = final equatorial helium ion density of flux tube (m^-3)
% hyden = northern hemisphere hydrogen density profiles (m^-3)
% heden = northern hemisphere helium density profiles (m^-3)
% oxden = northern hemisphere oxygen density profiles (m^-3)
% hydenc = southern hemisphere hydrogen density profiles (m^-3)
% hedenc = southern hemisphere helium density profiles (m^-3)
% oxdenc = southern hemisphere oxygen density profiles (m^-3)

%Declare required global variables
global LON350 LON350C LAT350 LAT350C LATMATM DISMAT DFLM LSHV MLTV
global EDEN350 HMF2 NMF2 BETA ETA DELTA ZETA EDEN350C HMF2C NMF2C BETAC ETAC DELTAC ZETAC
global OXNT350 OXNT350C HYNT350 HYNT350C HENT350 HENT350C N2NT350 N2NT350C
global TMPN350 TMPN350C TMPNEX TMPNEXC TMPNSP TMPNSPC

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)
rats=0.1;                               %diffusive/chemical loss ratio at which to start upwards diffusive profile
idhf2=50;                               %height above HmF2 at which use diffusive equilibrium for O+ (km)

%Select out required field line co-ordinates 
latmatm=LATMATM(poi,:);
dismat=DISMAT(poi,:);
dflm=DFLM(poi,:);

%Estimate the equatorial H+ and He+ densities by assuming they are equal to the respective mean densities
hyeq=hytot./fltpr(LSHV(poi));           %H+ equatorial density
if hes==1
 heeq=hetot./fltpr(LSHV(poi));          %He+ equatorial density
else
 heeq=zeros(size(hyeq));
end

%Calculate required ionosphere electron density values for the two bases of the field line.
[eden,edenc,hmf2p,hmf2pc]=iriap9(timi,dismat,latmatm,EDEN350,HMF2,NMF2,BETA,ETA,DELTA,ZETA,EDEN350C,HMF2C,NMF2C,BETAC,ETAC,DELTAC,ZETAC);

%Calculate the required neutral parameters. 
[oxnt,hynt,hent,n2nt,tmpn,oxntc,hyntc,hentc,n2ntc,tmpnc]=modmsis4(lsh,lshcp,timi,dismat,latmatm,OXNT350,HYNT350,HENT350,N2NT350,TMPN350,TMPNEX,TMPNSP,OXNT350C,HYNT350C,HENT350C,N2NT350C,TMPN350C,TMPNEXC,TMPNSPC,hes);

%Calculate the electron and ion temperatures.
[tmpe,tmpi,tmpec,tmpic,tdiur,tdiurc]=temppaw34(poi,t0d,t0n,g0d,g0n,t0dc,t0nc,g0dc,g0nc,bh,dhr,lsh,LON350(poi,timi),LAT350(poi,timi),LON350C(poi,timi),LAT350C(poi,timi),dismat,latmatm,tmpn,tmpnc,timut,dno,year);
tmpt=tmpi+tmpe;                                             %total temperature
tmptc=tmpic+tmpec;                                          %total conjugate temperature

%Determine the gravity component along the field at each point of interest.
g0=gravmod(dismat,latmatm);

%Set the oxygen density equal to the IRI2001 electron density values. It is assumed
%that the oxygen density is equal to the electron density for all points below
%HmF2+idhf2 km (using the HmF2 value at the base of the field line). For points 
%above this diffusive equilibrium will be used (see plsd8 below).
oxden=eden;                           %oxygen density at points of interest
oxdenc=edenc;                         %conjugate oxygen density at points of interest
sv=ones(1,size(oxden,2));

%Use the diffusive equilibrium model to estimate the intial O+ density profile.
[oxden,oxdenc]=plsd1(oxden,oxdenc,hyeq*sv,hyeq*sv,heeq*sv,heeq*sv,dflm,dismat,latmatm,g0,tmpi,tmpe,tmpt,tmpic,tmpec,tmptc,(re+hmf2p+idhf2),(re+hmf2pc+idhf2));

%Use the diffusive equilibrium model to estimate He+ density profiles. 
if hes==1
 [heden,hedenc]=plsd3(LSHV(poi),oxden,oxdenc,hyeq*sv,hyeq*sv,heeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc,tdiur*sv,tdiurc*sv,f107*ones(size(oxden,1),1));
else
 heden=zeros(size(oxden));
 hedenc=heden;
end

%Use the diffusive equilibrium model to estimate H+ density profiles. 
[hyden,hydenc]=plsd2(oxden,oxdenc,heden,hedenc,hyeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc);

%Correct to first order the H+ and He+ profiles so that the total tube content for the calculated profile agrees with orignal value. 
[hyden,hydenc,hyeq,heden,hedenc,heeq,csa]=tbcck(hyden,hydenc,heden,hedenc,latmatm,dismat,dflm,hytot,hyeq,hetot,heeq,hes);

%Recalculate the oxygen density profile.
[oxden,oxdenc,shoxm,shoxmc]=plsd1(oxden,oxdenc,hyden,hydenc,heden,hedenc,dflm,dismat,latmatm,g0,tmpi,tmpe,tmpt,tmpic,tmpec,tmptc,(re+hmf2p+idhf2),(re+hmf2pc+idhf2));

%Recalculate the helium density.
if hes==1
 [heden,hedenc,shhem,shhemc]=plsd3(LSHV(poi),oxden,oxdenc,hyden,hydenc,heeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc,tdiur*sv,tdiurc*sv,f107*ones(size(oxden,1),1));
end

%Recalculate the hydrogen density.
[hyden,hydenc,shhym,shhymc]=plsd2(oxden,oxdenc,heden,hedenc,hyeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc);

%Modifiy the H+ and conjugate H+ if required to include the effects of chemical equilibrium, and determine the downwards H+ flux.
[hyden,hyflux,hydiso,hydise]=plsd4(hyden,oxden,shhym,shoxm,dismat,dflm,tmpi,tmpe,tmpt,tmpn,oxnt,hynt,g0,r0,rats);
[hydenc,hyfluxc,hydisoc,hydisec]=plsd4(hydenc,oxdenc,shhymc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,oxntc,hyntc,g0,r0,rats);

%Modifiy the He+ and conjugate He+ if required to include the effects of chemical equilibrium, and determine the downwards H+ flux.
if hes==1
 [heden,heflux,hediso,hedise]=plsd6(heden,oxden,shhem,shoxm,dismat,dflm,tmpi,tmpe,tmpt,tmpn,n2nt,hent,g0,tdiur*sv,f107,f107a,r0,rats);
 [hedenc,hefluxc,hedisoc,hedisec]=plsd6(hedenc,oxdenc,shhemc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,n2ntc,hentc,g0,tdiurc*sv,f107,f107a,r0,rats);
end

%Calculate the nett upwards H+ flux and the diffusive downwards flux
[hyplsf,hyfluxt]=plsd5(LSHV(poi),latmatm,dflm,tmpi,tmpn,hyden,oxden,oxnt,hynt,shoxm,csa,hyflux,hydiso,hydise,kp,f107);
[hyplsfc,hyfluxtc]=plsd5(LSHV(poi),latmatm,dflm,tmpic,tmpnc,hydenc,oxdenc,oxntc,hyntc,shoxmc,csa,hyfluxc,hydisoc,hydisec,kp,f107);

%Determine the geographic longitude and latitude of the points on interest
lonmatm=timml2lon(timut,MLTV(poi),dno);
[lonmat,latmat]=mag2geo3(lonmatm*sv,latmatm);          %geographic co-ordinates of points of interest
[lonmatc,latmatc]=mag2geo3(lonmatm*sv,-latmatm);       %conjugate geographic co-ordinates of points of interest

%Calculate the nett upwards He+ flux and the diffusive downwards flux
[day,month]=daynum2(dno,year);
if hes==1
 [heplsf,hefluxt]=plsd7(LSHV(poi),dismat,latmatm,dflm,csa,lonmat,latmat,timut,heden,n2nt,hent,g0,tmpn,heflux,hediso,hedise,kp,f107,f107a,year,month,day);
 [heplsfc,hefluxtc]=plsd7(LSHV(poi),dismat,latmatm,dflm,csa,lonmatc,latmatc,timut,hedenc,n2ntc,hentc,g0,tmpnc,hefluxc,hedisoc,hedisec,kp,f107,f107a,year,month,day);
end

%% Determine new flux tube content
hytotn=hytot+3600*tint.*((hyplsf+hyplsfc)+(hyfluxt+hyfluxtc));
if hes==1
 hetotn=hetot+3600*tint.*((heplsf+heplsfc)+(hefluxt+hefluxtc));
end

%Check that the tube contents have not fallen below the 'background' value. This
%may occur if initial tube contents are unrealistically large and time interval
%choosen is too large.
[hyeqbg,heeqbg,hytotbg,hetotbg]=bgden(LSHV(poi),kp);
poibg=hytotn<hytotbg;
hytotn(poibg)=hytotbg(poibg);
if hes==1
 poibg2=hetotn<hetotbg;
 hetotn(poibg2)=hetotbg(poibg2);
end

%Equatorial densities
dvl=size(hyden,2);
hyeqf=hyden(:,dvl);
oxeqf=oxden(:,dvl);
if hes==1
 heeqf=heden(:,dvl);
else
 hetotn=zeros(size(hytotn));
 heeqf=zeros(size(hyeqf));
end
