%% plasdyn1.m is the improved version of the earlier plasdyn9.m,
%% which now includes He+ as well as O+ and H+. Further, it can 
%% also do the job of the earlier plasprof7.m.
%%
%% plasdyn1.m is self contained and should only be used when only 
%% a small number of field lines are required, as it uses the 
%% MSIS and IRI each time it is called.
%%
%% 5/3/2002 - He+ inclusion now optional. 
%%

function [hytotn,hetotn,hyden,oxden,heden,hydenc,oxdenc,hedenc,latmatm,dismat]=plasdyn1(lsh,hyinp,heinp,lonlsh,timmat,year,month,day,timi,nop,iricf,gms,iont,spp,hes,sms)

% lsh = L-shells of interest
% hyinp = total H+ content of flux tube or equatorial H+ density (m^-3)(set by inot)
% heinp = total He+ content of flux tube or equatorial He+ density (m^-3)(set by inot)
% lonlsh = magnetic longitude index of L-shells (rads) 
% timmat = univerisal time (fractional hours) (scalar)
% year = year (scalar)
% month = month (scalar)
% day = day (scalar)
% timi = time increment (fractional hours) (scalar)
% nop = number of points per field line (scalar)
% iricf = IRI2001 foF2 data set to be used (0 = CCIR, other = URSI) (scalar)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (scalar)
% iont = ion input type (0 = flux tube content, 1 = equatorial density) (scalar)
% spp = save parameters option (0 = no save, 1 = save) (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% hytotn = new total hydrogen ion content of flux tube (0 if inputs were equatorial densities)
% hetotn = new total helium ion content of flux tube (0 if inputs were equatorial densities)
% hyden = northern hemisphere hydrogen density profile (m^-3)
% oxden = northern hemisphere oxygen density profile (m^-3)
% heden = northern hemisphere helium density profile (m^-3)
% hydenc = southern hemisphere hydrogen density profile (m^-3)
% oxdenc = southern hemisphere oxygen density profile (m^-3)
% hedenc = southern hemisphere helium density profile (m^-3)
% latmatm = magnetic latitude of profile points (all positive)
% dismat = radial displacement of profile points (km)

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)

sz=size(hyinp);                         %size of input

%convert all input to column vectors
lsh=lsh(:);
hyinp=hyinp(:);
if hes==1
 heinp=heinp(:);
end
lonlsh=lonlsh(:);

lsh(lsh>15)=15;                         %reset terms with L-shell greater than 15 to 15
vlsh=fltpr(lsh);                        %volume of the L-shell (m^3)

%Set equatorial densities. If input was tube content then estimate the equatorial H+ and He+ 
%densities by assuming they are equal to the respective mean densities
if iont==0
 hyeq=hyinp./vlsh;
 hytot=hyinp;
 if hes==1
  heeq=heinp./vlsh;
  hetot=heinp;
 else
  heeq=zeros(size(hyeq));
 end
else
 hyeq=hyinp;
 if hes==1
  heeq=heinp;
 else
  heeq=zeros(size(hyeq));
 end
end

%Calculate the co-ordinate parameters of the field line
[lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm]=flpar(lsh,lonlsh,r0,nop);
sv=ones(1,size(lonmat,2));
sv2=ones(size(lonmat,1),1);
svm=ones(size(lonmat));

%Calculate required ionosphere electron density values for the two bases of the field line.
sv3=sv2*[1 1];
[eden350,nmf2,hmf2,beta,eta,delta,zeta]=modiri5([lonmat(:,1) lonmatc(:,1)],[latmat(:,1) latmatc(:,1)],timmat*sv3,day*sv3,month*sv3,year*sv3,r0*sv3,1000,iricf,gms);
[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 required solar parameters
dno=daynum(day,month,year);
[kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timmat,dno,year);

%Calculate the required neutral parameters.
svm2=[svm;svm];

apm=zeros([2*prod(sz) nop 7]);
for k=1:7
 apm(:,:,k)=apmat(k)*svm2;
end

[oxntt,hyntt,hentt,n2ntt,tmpnt]=modmsis1([lonmat;lonmatc],[latmat;latmatc],[dismat;dismat],timmat*svm2,day*svm2,month*svm2,year*svm2,f107a*svm2,f107p*svm2,apm,sms);
oxntt(oxntt<1)=1;
poi1=(1:size(lonmat,1))';
poi2=((size(lonmat,1)+1):(2*size(lonmat,1)))';
oxnt=oxntt(poi1,:);
hynt=hyntt(poi1,:);
tmpn=tmpnt(poi1,:);
oxntc=oxntt(poi2,:);
hyntc=hyntt(poi2,:);
tmpnc=tmpnt(poi2,:);
if hes==1
 hent=hentt(poi1,:);
 n2nt=n2ntt(poi1,:);
 hentc=hentt(poi2,:);
 n2ntc=n2ntt(poi2,:);
end

%Calculate the temperatures at the points of interest and conjugate points.
[tmpi,tmpe,tmpic,tmpec,tdiur,tdiurc]=temppaw2(lonmat,latmat,dismat,timmat*svm,year*svm,month*svm,day*svm,tmpn,tmpnc,kp*svm,f107p*svm);
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 IRI 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);


%If the input were tube contents then correct to first order the H+ and He+ profiles so that the total 
%tube content for the calculated profile agrees with orignal value. Otherwise, just determine the
%required cross-sectional areas.
if iont==0
 [hyden,hydenc,hyeq,heden,hedenc,heeq,csa]=tbcck(hyden,hydenc,heden,hedenc,latmatm,dismat,dflm,hytot,hyeq,hetot,heeq,hes);
else
 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./dismat).^3.*sqrt(1+3*sin(latmatm).^2);           %magnetic field of points of interest (T)
 csa=fls./bf;                                                %cross sectional area of point of interest (m^2)
end

%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));
else
 shhem=[];
 shhemc=[];
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,hedenb,herat]=plsd6(heden,oxden,shhem,shoxm,dismat,dflm,tmpi,tmpe,tmpt,tmpn,n2nt,hent,g0,tdiur,f107,f107a,r0,rats);
 [hedenc,hefluxc,hedisoc,hedisec,hedenbc,heratc]=plsd6(hedenc,oxdenc,shhemc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,n2ntc,hentc,g0,tdiurc,f107,f107a,r0,rats);
else
 hediso=[];
 hedise=[];
 hedenb=[];
 herat=[];
 hedisoc=[];
 hedisec=[];
 hedenbc=[];
 heratc=[];
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);
else
 heplsf=[];
 hefluxt=[];
 hepdfi=[];
 helsfi=[];
 heplsfc=[];
 hefluxtc=[];
 hepdfic=[];
 helsfic=[];
end

%Determine new flux tube contents if inputs were tube content - otherwise output 0.
if iont==0
 hytotn=hytot+3600*timi.*((hyplsf+hyplsfc)+(hyfluxt+hyfluxtc));
 hytotn=reshape(hytotn,sz);
 if hes==1
  hetotn=hetot+3600*timi.*((heplsf+heplsfc)+(hefluxt+hefluxtc));
  hetotn=reshape(hetotn,sz);
 else
  hetotn=zeros(size(hytotn));
 end
else
 hytotn=0;
 hetotn=0;
end

if spp==1 
 %Save required parameters to file 
 save plasp 
end 
