%% plasimod1.m is the GPID global plasmasphere-ionosphere model. It determines the total 
%% hydrogen ion content along each field line specified by the user. It also estimates
%% the equatorial hydrogen and oxygen ion densities along the field lines. The results
%% are saved to files for later access when required. 
%%
%% plasimod1.m is a function file, which can be called by plasimod.m if desired. 
%% plasimod1.m can either start a new simulation or continue an existing one.
%%
%% 3/8/2001 - now call IRI2001 on every day, since now includes daily storm variations.
%% 6/11/2001 - updated to include He+.
%% 25/1/2002 - removed plasmapause formation velocity.
%% 25/1/2002 - increase the number of MLT intervals from 48 to 96.
%% 5/3/2002 - remove different files name extension for IRI CCIR or URSI parameters (store iricf in arrays).
%% 5/3/2002 - He+ inclusion now optional. 
%% 9/3/2002 - now call NRLMSIS2000 version of MSIS.
%% 18/4/2002 - now assumes all contuation files end with x.
%%

function plasimod1(dno,year,dayt,iricf,gms,hes,ctr)

% dno = Starting day number of run (scalar)
% year = Starting year of run (scalar)
% dayt = Number of days to run (scalar)
% iricf = IRI coefficents to use (0 = CCIR, 1 = URSI)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
% ctr = continue existing run (0 = no, 1 = yes)

%Either load in previous array data (ctr==1) or determine new arrays.
if ctr==0
 %Fix the remaining parameters for now (change in later version?)
 nop=50;       %Number of points to break field line into.  
 notp=48;      %Number of magnetic local time (MLT) intervals to be used.
 tint=1;       %UT time iteration step size to be used (hours)

 sfs=1;        %Use Liemohn scaling factor with McIlwain E5D E-field model. (0 = original McIllwain, 1 = Liemohn) (scalar)                  
 sms=1;        %Determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar).

 re=6371;                                               %Earth radi (km)
 lshmx=9;                                               %maximum L-shell    
 lshcp=6.8;                                             %maximum L-shell for calculation of profiles 
 lsh=[(re+400)/re 1.1 1.15 1.2:0.1:6.8 7:0.2:lshmx];    %required L-shell vector  

 noti=24/tint;                                          %number of time increments to be used
 timmat=0:tint:(24-tint);                               %UT time vector
 mlt=(0:(24/notp):(24-(24/notp)));                      %magnetic local time vector
 
 [day,month]=daynum2(dno,year);                         %first day 
 
else
   
 %% Load in previous data file
 yrs=int2str(year);
 dns=int2str(dno);
 if dno<10
  dns=['00',dns];
 elseif dno>=10 & dno<100
  dns=['0',dns];
 end
 eval(['load d',yrs,'d',dns,'x'])
 
 %Determine required starting densities
 hytotj=hytot(:,end); 
 hyeqj=hyeq(:,end);
 if hes==1
  hetotj=hetot(:,end); 
  heeqj=heeq(:,end);
 else
  hetotj=zeros(size(hytotj));
  heeqj=zeros(size(hyeqj));
 end

 %First day to start simulation
 [dno,year]=daynum3(dno,year,1);
 [day,month]=daynum2(dno,year);
 
 re=6371;                                               %Earth radi (km)
 noti=24/tint;                                          %number of time increments to be used
 timmat=0:tint:(24-tint);                               %UT time vector
    
end   

%Calculate the required lookup tables for the equatorial motion of the flux tubes. 
dfs=exist('plsimeq1.mat');
if dfs==2
 load plsimeq1
else 
 [lsht,mltt,ratt,kpt]=hytgrid1(lsh,mlt,tint*3600,sfs);      %required arrays
 save plsimeq1 lsht mltt ratt kpt sfs                       %save for later use
end 

%%Set up the grid co-ordinates that the model will use.
%Declare global variables for co-ordinates.
global LATMATM DISMAT DFLM LSHV MLTV LON350 LON350C LAT350 LAT350C

[MLTV,LSHV]=meshgrid(mlt,lsh);               %required L-shell and MLT matrices
LSHV=LSHV(:);                                %convert LSHV to column vector
MLTV=MLTV(:);                                %convert MLTV to a column vector
lsh=lsh(:);                                  %convert lsh to a column vector
r0=re+350;                                   %base radial distance of field lines                                             
poip=LSHV<=lshcp;                            %field lines for which to complete full field line calculations
nof=length(LSHV);                            %total number of field lines to be modelled
nof2=sum(poip);                              %number of points to undertake full calculations
svy=ones(nof,1);

%Calculate the required co-ordinates.
[LATMATM,DISMAT,DFLM,LON350,LON350C,LAT350,LAT350C]=flpar2(poip,LSHV,MLTV,timmat,r0,nop,1,sms);

%Declare required global variables for IRI2001 and MSIS-90 models
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

%% Set the initial H+ and He+ tube contents and equatorial densities equal to the background value when Kp=0
%% if starting a new run - otherwise use values from previous run (already loaded above)
if ctr==0
 [hydenj,hedenj,hytotj,hetotj]=bgden(LSHV,zeros(nof,1));
end

%% Determine the initial base temperatures and gradients
[t0do,t0no,g0do,g0no,sl,bh,dhr]=temppaw31(lsh,DISMAT,LATMATM);

%% Now loop through the required number of days, of which each day consists
%% of 24/tint time steps of tint duration (hours) each
for i=1:dayt
 tic
 
 %Calculate i'th day look-up ionospheric parameter tables usign the IRI2001 model 
 [EDEN350,EDEN350C,HMF2,HMF2C,NMF2,NMF2C,BETA,BETAC,ETA,ETAC,DELTA,DELTAC,ZETA,ZETAC]=iriap10(dno,year,timmat,nof2,iricf,gms);
 
 %Calculate i'th day look-up neutral parameter tables using the MSIS model
 [OXNT350,OXNT350C,HYNT350,HYNT350C,HENT350,HENT350C,N2NT350,N2NT350C,TMPN350,TMPN350C,TMPNEX,TMPNEXC,TMPNSP,TMPNSPC]=modmsis3(dno,year,timmat,nof2,sms);
  
 %Result matrices
 hytot=zeros(nof,(noti+1));                  %total H+ tube content
 hyeq=zeros(nof,noti);                       %H+ equatorial density
 hetot=zeros(nof,(noti+1));                  %total He+ tube content
 heeq=zeros(nof,noti);                       %He+ equatorial density
 oxeq=hyeq;                                  %O+ equatorial density
 
 hytot(:,1)=hytotj;                          %record previous total H+ tube content density
 if hes==1
  hetot(:,1)=hetotj;                         %record previous total He+ tube content density
 end
  
 %%Loop through each time increment
 for j=1:noti
  timut=timmat(j);                           %i'th univerisal time
  disp(['   UT=',int2str(timut),'  Dno=',int2str(dno),'  Year=',int2str(year),'  (Day ',int2str(i),' of ',int2str(dayt),')'])  
  
  %Determine new parameters if start of new 3-hour period (ie 0000-0300, 0300-0600 etc)
  %Determine the seasonally correction and daily F10.7 for the base temperatures
  %and gradient for the given day number using the UT=0.
  if rem(timut,3)==0                                   
   [kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timut,dno,year);         %solar parameters     
   if timut==0
    [t0d,t0dc,t0n,t0nc,g0d,g0dc,g0n,g0nc]=temppaw32(t0do,t0no,g0do,g0no,dno,sl,f107);
   end
   %modify base temperatures and gradients due to Kp
   [t0d,t0dc,t0n,t0nc,g0d,g0dc,g0n,g0nc]=temppaw33(t0d,t0dc,t0n,t0nc,g0d,g0dc,g0n,g0nc,kp,sl);      
  end

  %Determine which field lines will have a new content calculated. This will
  %be the field lines inside L=lshcp. Points outside L=lshcp will retain their
  %current tube content. For points inside L=lshcp use plasdyn3.m to calculate 
  %their new field line content.  
  hydenj=-ones(nof,1);                  %j'th H+ equatorial densities
  hedenj=hydenj;                        %j'th He+ equatorial densities
  oxdenj=hydenj;                        %j'th O+ equatorial densities
    
  if any(poip)
   [hytotj(poip),hetotj(poip),oxdenj(poip),hydenj(poip),hedenj(poip)]=plasdyn3(poip,hytotj(poip),hetotj(poip),t0d,t0n,g0d,g0n,t0dc,t0nc,g0dc,g0nc,sl,bh,dhr,lsh,j,timut,dno,year,tint,lshcp,f107,f107a,kp,hes); 
  end    
    
  %Rotate the field lines through tint hours and then use GRIDDATA to determine
  %the tube content at the original grid location.
  hytotj2=reshape(hytotj,length(lsh),length(mlt));
  hytotj2=hytgrid2(hytotj2,lsh,mlt,kp,lsht,mltt,ratt,kpt,1,1,0);   
  hytotj=hytotj2(:);
  if hes==1
   hetotj2=reshape(hetotj,length(lsh),length(mlt));
   hetotj2=hytgrid2(hetotj2,lsh,mlt,kp,lsht,mltt,ratt,kpt,1,1,1);   
   hetotj=hetotj2(:); 
  end

  %Record j'th iteration results
  hytot(:,j+1)=hytotj;                  %record H+ tube content   
  hyeq(:,j)=hydenj;                     %record equatorial H+ densities  
  oxeq(:,j)=oxdenj;                     %record equatorial O+ densities  
  if hes==1
   hetot(:,j+1)=hetotj;                 %record He+ tube content   
   heeq(:,j)=hedenj;                    %record equatorial He+ densities  
  end
 end 
 
 %% Save the results
 yrs=int2str(year);
 dns=int2str(dno);
 if dno<10
  dns=['00',dns];
 elseif dno>=10 & dno<100
  dns=['0',dns];
 end
 
 if hes==1
  eval(['save d',yrs,'d',dns,'x hytot hetot mlt lsh oxeq hyeq heeq MLTV LSHV tint nop lshcp gms sfs hes iricf sms'])
 else
  eval(['save d',yrs,'d',dns,'x hytot mlt lsh oxeq hyeq MLTV LSHV tint nop lshcp gms sfs hes iricf sms'])
 end
 
 tt=toc;                                     %time taken to simulate day
 disp(['Time taken to simulate daynumber ',int2str(dno),' was ',int2str(tt),' seconds'])
 
 %% Determine the new day, month, year and day number 
 [dno,year]=daynum3(dno,year,1);             %new day number and year
 [day,month]=daynum2(dno,year);              %new day and month number
 
end
