%% plasires.m calculates the H+ tube content, and the O+ and H+
%% equatorial densities for the nearest 2 UT (above and below 
%% required timmat) using results produced by plasimod#. If 
%% only scalar values for each of the input are given (ie only 
%% one output required) then the output arrays are two dimensional 
%% with size lsh x mlt. If the input are vectors, then the output 
%% are three dimensional with size lsh x mlt x number of input.
%% 
%% NOTE - for multiple input the fastest results obtained with 
%%        results in chronological order if more than one result
%%        required for same day number.
%%
%% 9/10/00 - modify plasires.m so that it loads in an array called
%%           datm (see code below) 
%%           When plasires is called, it uses datm to see if the 
%%           required file is in the subdirectory GPIDres. If it
%%           is, it loads in the file from its sub directory. If
%%           the file is not listed in the datm, plasires assumes 
%%           it is in the same directory as plasires (which is how
%%           it use to work).
%%
%%           NOTE - any updates also need to be added to plasires2.m. 
%%
%% 10/11/2001 - update to include He+.
%% 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. 
%%

function [hytotil,hetotil,hyeqil,heeqil,oxeqil,hytotiu,hetotiu,hyeqiu,heeqiu,oxeqiu,timmatl,timmatu,lshi,mlti,lshcp,gms,iricf,hes]=plasires(timmat,dno,year)

% timmat = UT (scalar or vector)
% dno = day number (scalar or vector)
% year = year (scalar or vector)
%
% hytotil = lower H+ tube content output
% hetotil = lower He+ tube content output
% hyeqil = lower H+ equatorial density output
% heeqil = lower He+ equatorial density output
% oxeqil = lower O+ equatorial density output
% hytotiu = upper H+ tube content output
% hetotiu = upper He+ tube content output
% hyeqiu = upper H+ equatorial density output
% heeqiu = upper He+ equatorial density output
% oxeqiu = upper O+ equatorial density output
% timmatl = lower UT values vector
% timmatu = upper UT values vector
% lshi = L-shell vector
% mlti = MLT vector
% lshcp = maximum L-shell for calculation of profiles
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (scalar)
% iricf = IRI FoF2 data set used (0 = CCIR, 1 = URSI) (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)


%Known data sets previously calculated.
datm=plasiresdat;

%Load in plasimod data set for 1st required day
dnoj=dno(1);
yearj=year(1);
timmatj=timmat(1);

yrs=int2str(yearj);
dns=int2str(dnoj);
if dnoj<10
 dns=['00',dns];
elseif dnoj>=10 & dnoj<100
 dns=['0',dns];
end

%Determine if required data set is in GPIDres subdirectory. If it is,
%load file from directory, otherwise assume in current directory.
gpc=find(yearj==datm(:,1) & dnoj>=datm(:,3) & dnoj<=datm(:,4));

if isempty(gpc)
 eval(['load d',yrs,'d',dns])
else
 eval(['load GPIDres\y',int2str(datm(gpc,1)),'\run',int2str(datm(gpc,2)),'\d',yrs,'d',dns])
end

%% Determine required time indexes
timv=(0:tint:24);
timmx=(24/tint)+1;
timi1=sum(timmatj>=timv);
timi2=timi1+1;
timi2(timi2==(timmx+1))=timmx;

%% Select out required data sets and combine using weighting determined by UT
hytotil=hytot(:,timi1); 
hyeqil=hyeq(:,timi1);
oxeqil=oxeq(:,timi1);
hytotiu=hytot(:,timi2); 
hyeqiu=hyeq(:,timi2);
oxeqiu=oxeq(:,timi2);
timmatl=timv(timi1);
timmatu=timv(timi2);
if hes==1
 hetotil=hetot(:,timi1); 
 heeqil=heeq(:,timi1);
 hetotiu=hetot(:,timi2); 
 heeqiu=heeq(:,timi2);
end

%% Reshape to original size
ly=length(lsh);
lx=length(mlt);
hytotil=reshape(hytotil,ly,lx); 
hyeqil=reshape(hyeqil,ly,lx);
oxeqil=reshape(oxeqil,ly,lx);
hytotiu=reshape(hytotiu,ly,lx); 
hyeqiu=reshape(hyeqiu,ly,lx);
oxeqiu=reshape(oxeqiu,ly,lx);
if hes==1
 hetotil=reshape(hetotil,ly,lx); 
 heeqil=reshape(heeqil,ly,lx);
 hetotiu=reshape(hetotiu,ly,lx); 
 heeqiu=reshape(heeqiu,ly,lx);
end

%% Update look up data files to include a 24 MLT column using the 0 MLT values
hytotil=[hytotil hytotil(:,1)];
hyeqil=[hyeqil hyeqil(:,1)];
oxeqil=[oxeqil oxeqil(:,1)];
hytotiu=[hytotiu hytotiu(:,1)];
hyeqiu=[hyeqiu hyeqiu(:,1)];
oxeqiu=[oxeqiu oxeqiu(:,1)];
mlti=[mlt 24];
if hes==1
 hetotil=[hetotil hetotil(:,1)];
 heeqil=[heeqil heeqil(:,1)];
 hetotiu=[hetotiu hetotiu(:,1)];
 heeqiu=[heeqiu heeqiu(:,1)];
end

%% Update look up data files to include a L=(re+350)/re row. Assume it is 
%% equal to the lsh(1) values which is L=(re+400)/re.
re=6371;                           %Earth radi (km)
hytotil=[hytotil(1,:);hytotil];
hyeqil=[hyeqil(1,:);hyeqil];
oxeqil=[oxeqil(1,:);oxeqil];
hytotiu=[hytotiu(1,:);hytotiu];
hyeqiu=[hyeqiu(1,:);hyeqiu];
oxeqiu=[oxeqiu(1,:);oxeqiu];
lshi=[(re+350)/re;lsh];
if hes==1
 hetotil=[hetotil(1,:);hetotil];
 heeqil=[heeqil(1,:);heeqil];
 hetotiu=[hetotiu(1,:);hetotiu];
 heeqiu=[heeqiu(1,:);heeqiu];
end

%% If more than one set of results required, then loop through the remaining input.
nop=length(timmat);                %number of points
if nop>1
 hytotjl=zeros((ly+1),(lx+1),nop);      %total H+ content array
 hyeqjl=hytotjl;                        %equatorial H+ array
 oxeqjl=hytotjl;                        %equatorial O+ array
 hytotju=hytotjl;                       %total H+ content array
 hyeqju=hytotjl;                        %equatorial H+ array
 oxeqju=hytotjl;                        %equatorial O+ array
 if hes==1
  hetotjl=hytotjl;                      %equatorial He+ array
  heeqjl=hetotjl;                       %equatorial He+ array
  hetotju=hetotjl;                      %total He+ content array
  heeqju=hetotjl;                       %equatorial He+ array
 end
 
 hytotjl(:,:,1)=hytotil;                %record 1st results
 hyeqjl(:,:,1)=hyeqil;
 oxeqjl(:,:,1)=oxeqil;
 hytotju(:,:,1)=hytotiu;                %record 1st results
 hyeqju(:,:,1)=hyeqiu;
 oxeqju(:,:,1)=oxeqiu;
 if hes==1
  hetotjl(:,:,1)=hetotil;
  heeqjl(:,:,1)=heeqil;
  hetotju(:,:,1)=hetotiu;
  heeqju(:,:,1)=heeqiu;
 end
 
 %Rename 1st day and year loaded
 dnoj2=dnoj;
 yearj2=yearj;
 
 %Now through remaining input values
 for j=2:nop
    
  %% Load in plasimod data set for j'th required day or use previous day if same day
  dnoj=dno(j);
  yearj=year(j);
  timmatj=timmat(j);
  
  %If same day, use previously loaded data. Otherwise, load new data set
  if ~(dnoj==dnoj2 & yearj==yearj2)               
   yrs=int2str(yearj);
   dns=int2str(dnoj);
   if dnoj<10
    dns=['00',dns];
   elseif dnoj>=10 & dnoj<100
    dns=['0',dns];
   end

   %Determine if required data set is in GPIDres subdirectory. If it is,
   %load file from directory, otherwise assume in current directory.
   gpc=find(yearj==datm(:,1) & dnoj>=datm(:,3) & dnoj<=datm(:,4));
 
   if isempty(gpc)
    eval(['load d',yrs,'d',dns])
   else
    eval(['load GPIDres\y',int2str(datm(gpc,1)),'\run',int2str(datm(gpc,2)),'\d',yrs,'d',dns])
   end

  end
  
  %% Determine required time indexes
  timi1=sum(timmatj>=timv);
  timi2=timi1+1;
  timi2(timi2==(timmx+1))=timmx;

  %% Select out required data sets and combine using weighting determined by UT
  hytotil=hytot(:,timi1); 
  hyeqil=hyeq(:,timi1);
  oxeqil=oxeq(:,timi1);
  hytotiu=hytot(:,timi2); 
  hyeqiu=hyeq(:,timi2);
  oxeqiu=oxeq(:,timi2);
  timmatl(j)=timv(timi1);
  timmatu(j)=timv(timi2);
  if hes==1
   hetotil=hetot(:,timi1); 
   heeqil=heeq(:,timi1);
   hetotiu=hetot(:,timi2); 
   heeqiu=heeq(:,timi2);
  end
  
  %% Reshape to original size
  hytotil=reshape(hytotil,ly,lx); 
  hyeqil=reshape(hyeqil,ly,lx);
  oxeqil=reshape(oxeqil,ly,lx);
  hytotiu=reshape(hytotiu,ly,lx); 
  hyeqiu=reshape(hyeqiu,ly,lx);
  oxeqiu=reshape(oxeqiu,ly,lx);
  if hes==1
   hetotil=reshape(hetotil,ly,lx); 
   heeqil=reshape(heeqil,ly,lx);
   hetotiu=reshape(hetotiu,ly,lx); 
   heeqiu=reshape(heeqiu,ly,lx);
  end

  %% Update look up data files to include a 24 MLT column using the 0 MLT values
  hytotil=[hytotil hytotil(:,1)];
  hyeqil=[hyeqil hyeqil(:,1)];
  oxeqil=[oxeqil oxeqil(:,1)];
  hytotiu=[hytotiu hytotiu(:,1)];
  hyeqiu=[hyeqiu hyeqiu(:,1)];
  oxeqiu=[oxeqiu oxeqiu(:,1)];
  if hes==1
   hetotil=[hetotil hetotil(:,1)];
   heeqil=[heeqil heeqil(:,1)];
   hetotiu=[hetotiu hetotiu(:,1)];
   heeqiu=[heeqiu heeqiu(:,1)];
  end

  %% Update look up data files to include a L=(re+350)/re row. Assume it is 
  %% equal to the lsh(1) values which is L=(re+400)/re.
  hytotil=[hytotil(1,:);hytotil];
  hyeqil=[hyeqil(1,:);hyeqil];
  oxeqil=[oxeqil(1,:);oxeqil];
  hytotiu=[hytotiu(1,:);hytotiu];
  hyeqiu=[hyeqiu(1,:);hyeqiu];
  oxeqiu=[oxeqiu(1,:);oxeqiu];
  if hes==1
   hetotil=[hetotil(1,:);hetotil];
   heeqil=[heeqil(1,:);heeqil];
   hetotiu=[hetotiu(1,:);hetotiu];
   heeqiu=[heeqiu(1,:);heeqiu];
  end  

  %% Record j'th results
  hytotjl(:,:,j)=hytotil;         
  hyeqjl(:,:,j)=hyeqil;
  oxeqjl(:,:,j)=oxeqil;
  hytotju(:,:,j)=hytotiu;     
  hyeqju(:,:,j)=hyeqiu;
  oxeqju(:,:,j)=oxeqiu;
  if hes==1
   hetotjl(:,:,j)=hetotil;                 
   heeqjl(:,:,j)=heeqil;
   hetotju(:,:,j)=hetotiu;             
   heeqju(:,:,j)=heeqiu;
  end

  %% Update dnoj2 and yearj2
  dnoj2=dnoj;
  yearj2=yearj;
  
 end 
 
 %Rename output to required names
 hytotil=hytotjl;
 hyeqil=hyeqjl;
 oxeqil=oxeqjl;
 hytotiu=hytotju;
 hyeqiu=hyeqju;
 oxeqiu=oxeqju;
 if hes==1
  hetotil=hetotjl;
  heeqil=heeqjl;
  hetotiu=hetotju;
  heeqiu=heeqju;
 end 

end 

% Output empty He+ arrays for hes==0
if hes==0
 hetotil=[];
 heeqil=[];
 hetotiu=[];
 heeqiu=[];
end