%% plasires3.m calculates estimates of the H+ and He+ tube content, and the 
%% O+, H+ and He+ equatorial densities for the required points of interest. 
%% plasires3.m uses the results from plasimod(#) to determine the values
%% at the nearest UT above and below the required UT of the point, and 
%% then uses orbpath.m to determine the point of interest location at
%% these UT, where interpolation is used to determine their values. Final
%% the value at the required UT is obtained by taking the weighted mean 
%% of the two values based on their UTs and the point of interest UT.
%% 
%% NOTE - if all calculations are to undertaken at the same time for multiple
%%        points of interest, then the user should input scalar values for
%%        timmat, dno and year. This will give faster calculations.
%%
%% 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 [hytot,hetot,hyeq,heeq,oxeq,iricf,hes]=plasires3(timmat,dno,year,mlt,lsh)

% timmat = UT of points of interest (array or scalar)
% dno = day number of points of interest (array or scalar)
% year = year of points of interest (array or scalar)
% mlt = MLT of points of interest (array)
% lsh = L-shell of points of interest (array) 
%
% hytot = H+ tube content output (array)
% hetot = He+ tube content output (array)
% hyeq = H+ equatorial density output (array)
% heeq = He+ equatorial density output (array)
% oxeq = O+ equatorial density output (array)
% iricf = IRI2001 FoF2 data set used (0 = CCIR, other = URSI) (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)

sz=size(lsh);

%Convert all input to column vectors
timmat=timmat(:);
mlt=mlt(:);
lsh=lsh(:);
dno=dno(:);
year=year(:);

nol=length(timmat);
nop=length(lsh);

%Call plasires to determine values for adjacent UT of points of interest
[hytot1,hetot1,hyeq1,heeq1,oxeq1,hytot2,hetot2,hyeq2,heeq2,oxeq2,timmat1,timmat2,lshv,mltv,lshcp,gms,iricf,hes]=plasires(timmat,dno,year);

lshmx=max(lshv);

%Determine required Kp
kp=solarval(timmat,dno,year);

%%Loop through each point of interest. Determine it location at the
%%lower and upper UTs, use interpolation to determined required values
%%and value at required UT.

%First consider case of multiple points of interest at a single time
if nol==1
 dts1=3600*(timmat-timmat1);            %time to lower UT (secs)
 dts1=50*round(dts1/50);                %round to nearest multiple of 50
 dts2=3600*(timmat2-timmat);            %time to upper UT (secs)
 dts2=50*round(dts2/50);                %round to nearest multiple of 50
  
 %% Determine the H+ and He+ equatorial densities for each L-shell, including points outside
 %% the plasmapause. Therefore need to replace the -1 terms in the look up data tables
 %% with the density calculated from tube content divided by tube volume. Also calculate
 %% the O+ equatorial density, but as these values are only required to lshcp, don't
 %% need to replace -1 terms before doing the lookup.
 poih=hyeq1==-1;                        %terms to be replaced
 [poiy,poix]=find(poih);                %index of terms to be replaced
 volt=fltpr(lshv(poiy));                %volume of flux tubes to be replace
 hyeq1(poih)=hytot1(poih)./volt;        %lower H+ equatorial densities at -1 points
 hyeq2(poih)=hytot2(poih)./volt;        %upper H+ equatorial densities at -1 points
 if hes==1
  heeq1(poih)=hetot1(poih)./volt;       %lower He+ equatorial densities at -1 points
  heeq2(poih)=hetot2(poih)./volt;       %upper He+ equatorial densities at -1 points
 end

 %If UT of interest close enough to plasimod.m results, use values directly
 if dts1==0     
  hyeqpi=interp2(mltv,lshv,log(abs(hyeq1)),mlt,lsh);        %plasmasphere equatorial H+ densities
  hyeqpi=exp(hyeqpi);
  oxeqpi=interp2(mltv,lshv,log(abs(oxeq1)),mlt,lsh);        %plasmasphere equatorial O+ densities
  oxeqpi=exp(oxeqpi);
  hytoti=interp2(mltv,lshv,log(abs(hytot1)),mlt,lsh);       %plasmasphere H+ tube content
  hytoti=exp(hytoti);
  if hes==1
   heeqpi=interp2(mltv,lshv,log(abs(heeq1)),mlt,lsh);       %plasmasphere equatorial He+ densities
   heeqpi=exp(heeqpi);
   hetoti=interp2(mltv,lshv,log(abs(hetot1)),mlt,lsh);      %plasmasphere He+ tube content
   hetoti=exp(hetoti);
  end

  %Record required results
  hytot=hytoti(:);                             
  hyeq=hyeqpi(:);
  oxeq=oxeqpi(:);
  if hes==1
   hetot=hetoti(:);                             
   heeq=heeqpi(:);
  end
  
 elseif dts2==0 
  hyeqpi=interp2(mltv,lshv,log(abs(hyeq2)),mlt,lsh);        %plasmasphere equatorial H+ densities
  hyeqpi=exp(hyeqpi);
  oxeqpi=interp2(mltv,lshv,log(abs(oxeq2)),mlt,lsh);        %plasmasphere equatorial O+ densities
  oxeqpi=exp(oxeqpi);
  hytoti=interp2(mltv,lshv,log(abs(hytot2)),mlt,lsh);       %plasmasphere H+ tube content
  hytoti=exp(hytoti);
  if hes==1
   heeqpi=interp2(mltv,lshv,log(abs(heeq2)),mlt,lsh);       %plasmasphere equatorial He+ densities
   heeqpi=exp(heeqpi);
   hetoti=interp2(mltv,lshv,log(abs(hetot2)),mlt,lsh);      %plasmasphere He+ tube content
   hetoti=exp(hetoti);
  end  

  %Record required results
  hytot=hytoti(:);                             
  hyeq=hyeqpi(:);
  oxeq=oxeqpi(:);
  if hes==1
   hetot=hetoti(:);                             
   heeq=heeqpi(:);
  end
 
 %% If points of interest time is not close enough to the upper or lower UT, use
 %% orbpath.m to determine the points of interest locations at the upper and lower
 %% UT, use interpolation to determine densities and tube content at these values
 %% and then use weighting based on timmat to determine required values.
 else 
  [lsh1i,mlt1i]=orbpath(lsh,mlt,kp,dts1,0,1,1);             %location of point of interest in lower UT
  [lsh2i,mlt2i]=orbpath(lsh,mlt,kp,dts2,0,0,1);             %location of point of interest in upper UT
  
  %Check L-shell values have not moved outside allowed values
  lsh1i(lsh1i>lshmx)=lshmx;
  lsh2i(lsh2i>lshmx)=lshmx;
  
  %Interpolate lower values
  hyeqp1i=interp2(mltv,lshv,log(abs(hyeq1)),mlt1i,lsh1i);        %plasmasphere equatorial H+ densities
  hyeqp1i=exp(hyeqp1i);
  oxeqp1i=interp2(mltv,lshv,log(abs(oxeq1)),mlt1i,lsh1i);        %plasmasphere equatorial O+ densities
  oxeqp1i=exp(oxeqp1i);
  hytot1i=interp2(mltv,lshv,log(abs(hytot1)),mlt1i,lsh1i);       %plasmasphere H+ tube content
  hytot1i=exp(hytot1i);
  if hes==1
   heeqp1i=interp2(mltv,lshv,log(abs(heeq1)),mlt1i,lsh1i);       %plasmasphere equatorial He+ densities
   heeqp1i=exp(heeqp1i);
   hetot1i=interp2(mltv,lshv,log(abs(hetot1)),mlt1i,lsh1i);      %plasmasphere He+ tube content
   hetot1i=exp(hetot1i);
  end

  %Interpolate upper values
  hyeqp2i=interp2(mltv,lshv,log(abs(hyeq2)),mlt2i,lsh2i);        %plasmasphere equatorial H+ densities
  hyeqp2i=exp(hyeqp2i);
  oxeqp2i=interp2(mltv,lshv,log(abs(oxeq2)),mlt2i,lsh2i);        %plasmasphere equatorial O+ densities
  oxeqp2i=exp(oxeqp2i);
  hytot2i=interp2(mltv,lshv,log(abs(hytot2)),mlt2i,lsh2i);       %plasmasphere H+ tube content
  hytot2i=exp(hytot2i);
  if hes==1
   heeqp2i=interp2(mltv,lshv,log(abs(heeq2)),mlt2i,lsh2i);       %plasmasphere equatorial He+ densities
   heeqp2i=exp(heeqp2i);
   hetot2i=interp2(mltv,lshv,log(abs(hetot2)),mlt2i,lsh2i);      %plasmasphere He+ tube content
   hetot2i=exp(hetot2i);
  end

  %Determine required values as weighted mean of upper and lower UT values
  wt=(timmat-timmat1)/(timmat2-timmat1);
  hytot=(1-wt)*hytot1i(:)+wt*hytot2i(:);
  hyeq=(1-wt)*hyeqp1i(:)+wt*hyeqp2i(:);
  oxeq=(1-wt)*oxeqp1i(:)+wt*oxeqp2i(:);
  if hes==1
   hetot=(1-wt)*hetot1i(:)+wt*hetot2i(:);
   heeq=(1-wt)*heeqp1i(:)+wt*heeqp2i(:);
  end
 end
 
%% If separate time calculations are required, loop through each each point
%% of interest and calculate separately the required densities and tube content. 
else  
 %Results arrays
 hytot=zeros(nop,1);
 hyeq=hytot;
 oxeq=hytot;
 if hes==1
  hetot=hytot;                      
  heeq=hytot;
 end
 
 %Need to loop through each point of interest separately
 for j=1:nop
  dts1=3600*(timmat(j)-timmat1(j));          %time to lower UT (secs)
  dts1=50*round(dts1/50);                    %round to nearest multiple of 50
  dts2=3600*(timmat2(j)-timmat(j));          %time to upper UT (secs)
  dts2=50*round(dts2/50);                    %round to nearest multiple of 50
  
  %% Determine the H+ and He+ equatorial densities for each L-shell, including points outside
  %% the plasmapause. Therefore need to replace the -1 terms in the look up data tables
  %% with the density calculated from tube content divided by tube volume. Also calculate
  %% the O+ equatorial density, but as these values are only required to lshcp, don't
  %% need to replace -1 terms before doing the lookup.
  hyeq1j=squeeze(hyeq1(:,:,j));         %j'th lower H+ densities
  oxeq1j=squeeze(oxeq1(:,:,j));         %j'th lower O+ densities
  hytot1j=squeeze(hytot1(:,:,j));       %j'th lower H+ tube content
  if hes==1
   heeq1j=squeeze(heeq1(:,:,j));        %j'th lower He+ densities
   hetot1j=squeeze(hetot1(:,:,j));      %j'th lower He+ tube content
  end 

  hyeq2j=squeeze(hyeq2(:,:,j));         %j'th upper H+ densities
  oxeq2j=squeeze(oxeq2(:,:,j));         %j'th upper O+ densities
  hytot2j=squeeze(hytot2(:,:,j));       %j'th upper H+ tube content
  if hes==1
   heeq2j=squeeze(heeq2(:,:,j));        %j'th upper He+ densities
   hetot2j=squeeze(hetot2(:,:,j));      %j'th upper He+ tube content
  end
  
  poih=hyeq1j==-1;                      %terms to be replaced
  [poiy,poix]=find(poih);               %index of terms to be replaced
  volt=fltpr(lshv(poiy));               %volume of flux tubes to be replace
  hyeq1j(poih)=hytot1j(poih)./volt;     %lower H+ equatorial densities at -1 points
  hyeq2j(poih)=hytot2j(poih)./volt;     %upper H+ equatorial densities at -1 points
  if hes==1
   heeq1j(poih)=hetot1j(poih)./volt;    %lower He+ equatorial densities at -1 points
   heeq2j(poih)=hetot2j(poih)./volt;    %upper He+ equatorial densities at -1 points
  end

  %If UT of interest close enough to plasimod.m results, use values directly
  if dts1==0     
   hyeqpi=interp2(mltv,lshv,log(abs(hyeq1j)),mlt(j),lsh(j));     %plasmasphere equatorial H+ densities
   hyeqpi=exp(hyeqpi);
   oxeqpi=interp2(mltv,lshv,log(abs(oxeq1j)),mlt(j),lsh(j));     %plasmasphere equatorial O+ densities
   oxeqpi=exp(oxeqpi);
   hytoti=interp2(mltv,lshv,log(abs(hytot1j)),mlt(j),lsh(j));    %plasmasphere H+ tube content
   hytoti=exp(hytoti);
   if hes==1
    heeqpi=interp2(mltv,lshv,log(abs(heeq1j)),mlt(j),lsh(j));    %plasmasphere equatorial He+ densities
    heeqpi=exp(heeqpi);
    hetoti=interp2(mltv,lshv,log(abs(hetot1j)),mlt(j),lsh(j));   %plasmasphere He+ tube content
    hetoti=exp(hetoti);
   end

   %Record required results
   hytot(j)=hytoti;                             
   hyeq(j)=hyeqpi;
   oxeq(j)=oxeqpi;
   if hes==1
    hetot(j)=hetoti;                             
    heeq(j)=heeqpi;
   end
  
  elseif dts2==0 
   hyeqpi=interp2(mltv,lshv,log(abs(hyeq2j)),mlt(j),lsh(j));     %plasmasphere equatorial H+ densities
   hyeqpi=exp(hyeqpi);
   oxeqpi=interp2(mltv,lshv,log(abs(oxeq2j)),mlt(j),lsh(j));     %plasmasphere equatorial O+ densities
   oxeqpi=exp(oxeqpi);
   hytoti=interp2(mltv,lshv,log(abs(hytot2j)),mlt(j),lsh(j));    %plasmasphere H+ tube content
   hytoti=exp(hytoti);
   if hes==1
    heeqpi=interp2(mltv,lshv,log(abs(heeq2j)),mlt(j),lsh(j));    %plasmasphere equatorial He+ densities
    heeqpi=exp(heeqpi);
    hetoti=interp2(mltv,lshv,log(abs(hetot2j)),mlt(j),lsh(j));   %plasmasphere He+ tube content
    hetoti=exp(hetoti);
   end 
  
   %Record required results
   hytot(j)=hytoti;                             
   hyeq(j)=hyeqpi;
   oxeq(j)=oxeqpi;
   if hes==1
    hetot(j)=hetoti;                             
    heeq(j)=heeqpi;
   end
  
  %% If points of interest time is not close enough to the upper or lower UT, use
  %% orbpath.m to determine the points of interest locations at the upper and lower
  %% UT, use interpolation to determine densities and tube content at these values
  %% and then use weighting based on timmat to determine required values.
  else 
   [lsh1i,mlt1i]=orbpath(lsh(j),mlt(j),kp(j),dts1,0,1,1);        %location of point of interest in lower UT
   [lsh2i,mlt2i]=orbpath(lsh(j),mlt(j),kp(j),dts2,0,0,1);        %location of point of interest in upper UT
   
   %Check L-shell values have not moved outside allowed values
   lsh1i(lsh1i>lshmx)=lshmx;
   lsh2i(lsh2i>lshmx)=lshmx;
    
   %Interpolate lower values
   hyeqp1i=interp2(mltv,lshv,log(abs(hyeq1j)),mlt1i,lsh1i);      %plasmasphere equatorial H+ densities
   hyeqp1i=exp(hyeqp1i);
   oxeqp1i=interp2(mltv,lshv,log(abs(oxeq1j)),mlt1i,lsh1i);      %plasmasphere equatorial O+ densities
   oxeqp1i=exp(oxeqp1i);
   hytot1i=interp2(mltv,lshv,log(abs(hytot1j)),mlt1i,lsh1i);     %plasmasphere H+ tube content
   hytot1i=exp(hytot1i);  
   if hes==1
    heeqp1i=interp2(mltv,lshv,log(abs(heeq1j)),mlt1i,lsh1i);     %plasmasphere equatorial He+ densities
    heeqp1i=exp(heeqp1i);
    hetot1i=interp2(mltv,lshv,log(abs(hetot1j)),mlt1i,lsh1i);    %plasmasphere He+ tube content
    hetot1i=exp(hetot1i);
   end

   %Interpolate upper values
   hyeqp2i=interp2(mltv,lshv,log(abs(hyeq2j)),mlt2i,lsh2i);      %plasmasphere equatorial H+ densities
   hyeqp2i=exp(hyeqp2i);
   oxeqp2i=interp2(mltv,lshv,log(abs(oxeq2j)),mlt2i,lsh2i);      %plasmasphere equatorial O+ densities
   oxeqp2i=exp(oxeqp2i);
   hytot2i=interp2(mltv,lshv,log(abs(hytot2j)),mlt2i,lsh2i);     %plasmasphere H+ tube content
   hytot2i=exp(hytot2i);
   if hes==1
    heeqp2i=interp2(mltv,lshv,log(abs(heeq2j)),mlt2i,lsh2i);     %plasmasphere equatorial He+ densities
    heeqp2i=exp(heeqp2i);
    hetot2i=interp2(mltv,lshv,log(abs(hetot2j)),mlt2i,lsh2i);    %plasmasphere He+ tube content
    hetot2i=exp(hetot2i);
   end

   %Determine required values as weighted mean of upper and lower UT values
   wt=(timmat(j)-timmat1(j))/(timmat2(j)-timmat1(j));
   hytot(j)=(1-wt)*hytot1i+wt*hytot2i;
   hyeq(j)=(1-wt)*hyeqp1i+wt*hyeqp2i;
   oxeq(j)=(1-wt)*oxeqp1i+wt*oxeqp2i;
   if hes==1
    hetot(j)=(1-wt)*hetot1i+wt*hetot2i;
    heeq(j)=(1-wt)*heeqp1i+wt*heeqp2i;
   end
  end 
 end
end

%Reset any H+ and He+ densities below the background densities or equal to NaN to the background
[hyeqbg,heeqbg]=bgden(lsh,kp);
poihy=hyeq<hyeqbg | isnan(hyeq);
hyeq(poihy)=hyeqbg(poihy);
if hes==1
 poihe=heeq<heeqbg | isnan(heeq);
 heeq(poihe)=heeqbg(poihe);
end

%Reset any O+ densities below 1 to 1
oxeq(oxeq<1)=1;

%Reset any points greater than lshmx to background values, and O+ to 1.
poimx=lsh>lshmx;
hyeq(poimx)=hyeqbg(poimx);
hytot(poimx)=hyeqbg(poimx).*fltpr(lsh(poimx));
oxeq(poimx)=1;
if hes==1
 heeq(poimx)=heeqbg(poimx);
 hetot(poimx)=heeqbg(poimx).*fltpr(lsh(poimx));
end

% Reshape results to input size
hytot=reshape(hytot,sz);
hyeq=reshape(hyeq,sz);
oxeq=reshape(oxeq,sz);
if hes==1
 hetot=reshape(hetot,sz);
 heeq=reshape(heeq,sz);
end

% Output empty He+ arrays for hes==0
if hes==0
 hetot=[];
 heeq=[];
end