%% plasdyn2.m calculate the H+ and He+ refilling of a flux tube.
%%
%% For every 10 days in calls the IRI2001 and MSIS86 once only, storing
%% all required data in memory to help speed up excution. It only
%% returns the tube content and equatorial H+ and O+ densities.
%% Note that this approach assumes that the flux tube of interest is
%% co-rotating with the Earth and is not being emptied by magnetic
%% storms.
%%
%% 5/3/2002 - He+ inclusion now optional. 
%%

function [hytot,hetot,hyeq,heeq,oxeq,timm,dnom,yearm]=plasdyn2(lsh,hytoti,hetoti,lonlsh,nod,timi,year,month,day,timmat,nop,hynsf,hensf,iricf,gms,spp,hes,sms)

% lsh = L-shells of interest (scalar)
% hytoti = initial total H+ content of flux tube (scalar)
% hetoti = initial total He+ content of flux tube (scalar)
% lonlsh = magnetic longitude index of L-shells (rads) (scalar) 
% nod = number of days to run (scalar)
% timi = time increment (fractional hours) (scalar)
% year = starting year (scalar)
% month = starting month (scalar)
% day = starting day (scalar)
% timmat = starting UT (scalar)
% nop = number of points per field line (scalar)
% hynsf = neutral hydrogen scaling factor, all MSIS H densities multiplied by this number. 1 = no change (scalar)
%         NOTE - if set to minus value, then southern hemisphere value is taken as been 
%                equal to nothern value times abs(hynsf). 
% hensf = neutral helium scaling factor, all MSIS He densities multiplied by this number. 1 = no change (scalar)
%         NOTE - if set to minus value, then southern hemisphere value is taken as been 
%                equal to nothern value times abs(hensf). 
% iricf = IRI2001 FoF2 data set to be used (0 = CCIR, other = URSI) (scalar)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (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)
%
% hytot = total H+ content of flux tube 
% hetot = total He+ content of flux tube 
% hyeq = equatorial H+ density of flux tube (m^-3)
% heeq = equatorial He+ density of flux tube (m^-3)
% oxeq = equatorial O+ density of flux tube (m^-3)
% timm = UT time vector used
% dnom = day number vector used
% yearm = year vector used

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

noid=24/timi;                           %number of iterations per day
nodi=10;                                %number of days per iterative loop
noii=nodi*noid;                         %number of iterations to conduct per loop
noit=nod*noid;                          %total number of iterations required
noil=ceil(noit/noii);                   %total number of loops required

%Set up results arrays
hytot=zeros(noit,1);                    %total H+ content of flux tube
hetot=hytot;                            %total He+ content of flux tube
hyeq=hytot;                             %equatorial H+ density
heeq=hytot;                             %equatorial He+ density
oxeq=hytot;                             %equatorial O+ density
timm=hytot;                             %UT times used
dnom=hytot;                             %day numbers used
yearm=hytot;                            %year used

hyplsf=hytot;                           %nett upwards H+ production
hyplsfc=hytot;                          %nett conjugate H+ upwards production
hyflux=hytot;                           %downwards H+ flux
hyfluxc=hytot;                          %downwards H+ flux
hyplsfmx=hytot;                         %maximum upwards Richards and Torr H+ flux
hyplsfmxc=hytot;                        %maximum upwards Richards and Torr conjugate H+ flux
hychem=hytot;                           %nett upwards H+ production before Richards and Torr limit
hychemc=hytot;                          %nett conjugate H+ upwards production before Richards and Torr limit
hydise=hytot;                           %distance along field line where H+ chemical and diffusive loss equal
hydisec=hytot;                          %distance along field line where H+ conjugate chemical and diffusive loss equal
hytmpn=hytot;                           %neutral temperature at diso
hytmpnc=hytot;                          %conjugate neutral temperature at disoc
hyntd=hytot;                            %H density at diso
hyntdc=hytot;                           %conjugate H density at disoc
hyoxd=hytot;                            %O+ density at diso
hyoxdc=hytot;                           %conjugate O+ density at disoc
hyshoxd=hytot;                          %O+ scale height at diso
hyshoxdc=hytot;                         %conjugate O+ scale height at disoc

heplsf=hytot;                           %nett upwards He+ production
heplsfc=hytot;                          %nett conjugate He+ upwards production
heflux=hytot;                           %downwards He+ flux
hefluxc=hytot;                          %downwards He+ flux
heplsfo=hytot;                          %nett upwards He+ production
heplsfoc=hytot;                         %nett conjugate He+ upwards production
hedise=hytot;                           %distance along field line where He+ chemical and diffusive loss equal
hedisec=hytot;                          %distance along field line where He+ conjugate chemical and diffusive loss equal
hesh=hytot;                             %He+ scale height at grid point below hediso
heshc=hytot;                            %conjugate He+ scale height at grid point below hedisoc
oxsh=hytot;                             %O+ scale height at grid point below hediso
oxshc=hytot;                            %conjugate O+ scale height at grid point below hedisoc
n2sh=hytot;                             %N2 scale height at grid point below hediso
n2shc=hytot;                            %conjugate N2 scale height at grid point below hedisoc
hentsh=hytot;                           %He scale height at grid point below hediso
hentshc=hytot;                          %conjugate He scale height at grid point below hedisoc
hentd=hytot;                            %H density at grid point below hediso
hentdc=hytot;                           %conjugate H density at grid point below hedisoc
hehed=hytot;                            %H+ density at grid point below hediso
hehedc=hytot;                           %conjugate H+ density at grid point below hedisoc
hen2d=hytot;                            %N2 density at grid point below hediso
hen2dc=hytot;                           %conjugate N2 density at grid point below hedisoc
hedise=hytot;                           %distance along field line where He+ chemical and diffusive loss equal
hedisec=hytot;                          %distance along field line where He+ conjugate chemical and diffusive loss equal

%Set up required co-ordinates arrays
[lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm]=flpar(lsh,lonlsh,r0,nop);

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

%Loop through noil times, for each calculating the required IRI and MSIS parameters.
%Then within that, loop through noii times, calculating the evolution of the 
%flux tube H+ and He+ content.
dnoi=daynum(day,month,year);                 %initial day of i'th loop
yeari=year;                                  %initial year of i'th loop
[dnomx,yearmx]=daynum3(dnoi,yeari,nod);      %final daynumber and year

hyi=1;                                       %initial hytot index term

for i=1:noil
 
 %Determine the required time vectors for the i'th noil iteration
 if i~=noil
  nodi2=nodi;
  noii2=noii;
 else
  nodi2=rem(noit,noii)/noid;
  if nodi2==0
   nodi2=nodi;
  end
  noii2=nodi2*noid; 
 end
 
 tim0=timmat+(0:timi:((24*nodi2)-timi));     %total time vector
 tim0=tim0(:);                               %convert to column vector
 timv=rem(tim0,24);                          %required time vector for i'th loop
 dnoc=floor(tim0/24);                        %change in day number
 svt=ones(noii2,1);

 [dnov,yearv]=daynum3(dnoi*svt,yeari*svt,dnoc);   %required day number and year vector for i'th loop
 [dayv,monthv]=daynum2(dnov,yearv);               %required day and month vector for i'th loop
 
 nopt=length(timv);
 
 sv=ones(nopt,1);
 svm=ones(nopt,2);
 
 %Calculate the required solar parameters
 [kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timv,dnov,yearv);

 %Determine the required MSIS and IRI paramters
 [eden350,nmf2,hmf2,beta,eta,delta,zeta]=modiri5(sv*[lonmat(1) lonmatc(1)],sv*[latmat(1) latmatc(1)],timv*[1 1],dayv*[1 1],monthv*[1 1],yearv*[1 1],r0*svm,1000,iricf,gms);
 svx=ones(1,nop);
 apm=zeros([2*nopt nop 7]);
 for k=1:7
  apm(:,:,k)=[squeeze(apmat(:,:,k));squeeze(apmat(:,:,k))]*svx;
 end

 [oxntt,hyntt,hentt,n2ntt,tmpnt]=modmsis1([sv*lonmat;sv*lonmatc],[sv*latmat;sv*latmatc],[sv*dismat;sv*dismat],[timv;timv]*svx,[dayv;dayv]*svx,[monthv;monthv]*svx,[yearv;yearv]*svx,[f107a;f107a]*svx,[f107p;f107p]*svx,apm,sms);
 poi1=1:nopt;
 poi2=(nopt+1):(2*nopt);
 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
 
 %Scale the neutral H as required
 if hynsf>0
  hynt=hynsf*hynt;
  hyntc=hynsf*hyntc;
 else
  hyntc=abs(hynsf).*hynt;
 end
 
 %Scale the neutral He as required
 if hes==1
  if hensf>0
   hent=hensf*hent;
   hentc=hensf*hentc;
  else
   hentc=abs(hensf).*hent;
  end
 end
 
 %Loop through each separate time interval to determine the time evolution of the flux tube
 for k=1:length(timv) 
  if timv(k)<timi
   [dnov(k) yearv(k) dnomx yearmx]
  end 
  if hyi==1
   hytoto=hytoti;
   hyeqo=hytoti/fltpr(lsh);
   if hes==1
    hetoto=hetoti;
    heeqo=hetoti/fltpr(lsh);
   else
    hetoto=zeros(size(hytoto));
    heeqo=zeros(size(hyeqo));
   end
  else
   hytoto=hytot(hyi-1); 
   hyeqo=hyeq(hyi-1);
   if hes==1
    hetoto=hetot(hyi-1); 
    heeqo=heeq(hyi-1);
   end
  end  

  [hytotk,hetotk,hydenk,hedenk,oxdenk,hym,hymc,hem,hemc]=plasdyn21(hytoto,hetoto,hyeqo,heeqo,timv(k),dayv(k),monthv(k),yearv(k),lsh,lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm,eden350(k,:),nmf2(k,:),hmf2(k,:),beta(k,:),eta(k,:),delta(k,:),zeta(k,:),oxnt(k,:),hynt(k,:),hent(k,:),n2nt(k,:),tmpn(k,:),oxntc(k,:),hyntc(k,:),hentc(k,:),n2ntc(k,:),tmpnc(k,:),kp(k),f107(k),f107a(k),g0,timi,hes);
 
  hytot(hyi)=hytotk;
  hyeq(hyi)=hydenk;
  oxeq(hyi)=oxdenk;
  timm(hyi)=timv(k);
  dnom(hyi)=dnov(k);
  yearm(hyi)=yearv(k);
  if hes==1
   hetot(hyi)=hetotk;
   heeq(hyi)=hedenk;
  end

  hyplsf(hyi)=hym(1);
  hyplsfc(hyi)=hymc(1);
  hyflux(hyi)=hym(2);
  hyfluxc(hyi)=hymc(2);  
  hyplsfmx(hyi)=hym(3);
  hyplsfmxc(hyi)=hymc(3);
  hychem(hyi)=hym(4);
  hychemc(hyi)=hymc(4);
  hydise(hyi)=hym(5);
  hydisec(hyi)=hymc(5);
  hytmpn(hyi)=hym(6);
  hytmpnc(hyi)=hymc(6);
  hyntd(hyi)=hym(7);
  hyntdc(hyi)=hymc(7);
  hyoxd(hyi)=hym(8);
  hyoxdc(hyi)=hymc(8);
  hyshoxd(hyi)=hym(9);
  hyshoxdc(hyi)=hymc(9);

  heplsf(hyi)=hem(1);
  heplsfc(hyi)=hemc(1);
  heflux(hyi)=hem(2);
  hefluxc(hyi)=hemc(2);  
  hesh(hyi)=hem(3);
  heshc(hyi)=hemc(3);
  oxsh(hyi)=hem(4);
  oxshc(hyi)=hemc(4);
  n2sh(hyi)=hem(5);
  n2shc(hyi)=hemc(5);
  hentsh(hyi)=hem(6);
  hentshc(hyi)=hemc(6);
  hedise(hyi)=hem(7);
  hedisec(hyi)=hemc(7);
  hentd(hyi)=hem(8);
  hentdc(hyi)=hemc(8);
  hehed(hyi)=hem(9);
  hehedc(hyi)=hemc(9);
  hen2d(hyi)=hem(10);
  hen2dc(hyi)=hemc(10);
 
  hyi=hyi+1;
 end
 
 [dnoi,yeari]=daynum3(dnoi,yeari,nodi);
 
end

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



