% plasdyn21.m is designed to be called by plasdyn2.m, see it for input parameters description.
%

function [hytotn,hetotn,hyeqn,heeqn,oxeqn,hym,hymc,hem,hemc]=plasdyn21(hytot,hetot,hyeq,heeq,timmat,day,month,year,lsh,lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm,eden350,nmf2,hmf2,beta,eta,delta,zeta,oxnt,hynt,hent,n2nt,tmpn,oxntc,hyntc,hentc,n2ntc,tmpnc,kp,f107,f107a,g0,timi,hes)

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)

nopd=length(dismat);
sv=ones(1,nopd);

%Calculate required ionosphere electron density values for the two bases of the field line.
[eden,edenc]=plasdyn11(dismat,latmatm,eden350(1),hmf2(1),nmf2(1),beta(1),eta(1),delta(1),zeta(1),eden350(2),hmf2(2),nmf2(2),beta(2),eta(2),delta(2),zeta(2));

%Calculate the temperatures at the points of interest and conjugate points.
[tmpi,tmpe,tmpic,tmpec,tdiur,tdiurc]=temppaw2(lonmat,latmat,dismat,timmat*sv,year*sv,month*sv,day*sv,tmpn,tmpnc,kp*sv,f107*sv);
tmpt=tmpi+tmpe;                                             %total temperature
tmptc=tmpic+tmpec;                                          %total conjugate temperature

%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

%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+hmf2(:,1)+idhf2),(re+hmf2(:,2)+idhf2));

%Use the diffusive equilibrium model to estimate He+ density profiles. 
if hes==1
 [heden,hedenc]=plsd3(lsh,oxden,oxdenc,hyeq*sv,hyeq*sv,heeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc,tdiur,tdiurc,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+hmf2(:,1)+idhf2),(re+hmf2(:,2)+idhf2));

%Recalculate the helium density.
if hes==1
 [heden,hedenc,shhem,shhemc]=plsd3(lsh,oxden,oxdenc,hyden,hydenc,heeq,dflm,g0,tmpi,tmpe,tmpt,tmpn,tmpic,tmpec,tmptc,tmpnc,tdiur,tdiurc,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,f107,f107a,r0,rats);
 [hedenc,hefluxc,hedisoc,hedisec]=plsd6(hedenc,oxdenc,shhemc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,n2ntc,hentc,g0,tdiurc,f107,f107a,r0,rats);
end

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

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

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

%Equatorial densities
hyeqn=hyden(nopd);
if hes==1
 heeqn=heden(nopd);
end
oxeqn=oxden(nopd);

%Determine location of nearest grid point below diso
if hes==1
 pdis=sum(dflm<(hediso*sv));
 pdisc=sum(dflm<(hedisoc*sv));
end
pdis2=sum(dflm<(hydiso*sv));
pdis2c=sum(dflm<(hydisoc*sv));

%Output arrays
hym=zeros(1,9);
hymc=hym;
hym(1)=hyplsf;
hymc(1)=hyplsfc;
hym(2)=hyfluxt;
hymc(2)=hyfluxtc;
hym(3)=hyplsfmx;
hymc(3)=hyplsfmxc;
hym(4)=hychem;
hymc(4)=hychemc;
hym(5)=hydise;
hymc(5)=hydisec;
hym(6)=tmpn(pdis2);
hymc(6)=tmpnc(pdis2);
hym(7)=hynt(pdis2);
hymc(7)=hyntc(pdis2);
hym(8)=oxden(pdis2);
hymc(8)=oxdenc(pdis2);
hym(9)=shoxm(pdis2);
hymc(9)=shoxmc(pdis2);

hem=zeros(1,10);
hemc=hem;
if hes==1
 hem(1)=heplsf;
 hemc(1)=heplsfc;
 hem(2)=hefluxt;
 hemc(2)=hefluxtc;
 hem(3)=shhem(pdis);
 hemc(3)=shhemc(pdisc);
 hem(4)=shoxm(pdis);
 hemc(4)=shoxmc(pdisc);

 %Calculate neutral N2 and He scale heights pdis and pdisc points
 g0=9.8*(re./dismat(pdis)).^2;     
 g0c=9.8*(re./dismat(pdisc)).^2;                                  
                             
 hem(5)=1e-3*k*tmpn(pdis)./(28*amu*g0); 
 hemc(5)=1e-3*k*tmpnc(pdisc)./(28*amu*g0c); 

 hem(6)=7*hem(6); 
 hemc(6)=7*hemc(6); 

 hem(7)=hedise;
 hemc(7)=hedisec;
 hem(8)=hent(pdis);
 hemc(8)=hentc(pdisc);
 hem(9)=heden(pdis);
 hemc(9)=hedenc(pdisc);
 hem(10)=n2nt(pdis);
 hemc(10)=n2ntc(pdisc);
end

if hes==0
 hetotn=nan*ones(size(hytotn));
 heeqn=nan*ones(size(hyeqn));
end