%% plasieq2.m determines the electron densities for the specified magnetic
%% equatorial 'slice' at the given MLT, UT, day, month and year using data files
%% produced by plasimod(2).m. It is based on plasieq1.m, but returns arrays
%% in which each row corresponds to the required MLT 'radial ray' ie does not return
%% a single complete equatorial slice at one time. 
%%
%% NOTE - above lshcp for the data set, eden, oxden and hyden return NaN.
%%
%% NOTE - plasisc.m uses the equatorial values from the closest UT to the
%%        points of interest ie it does not interpolate between the nearest
%%        two times. 
%%

function [eden,oxden,hyden,heden,hytot,hetot,lshv]=plasieq2(mlt,timmat,dno,year)

% mlt = magnetic local time of radial path (hours)
% timmat = univerisal time (fractional hours)
% dno = day number
% year = year
%
% eden = electron density (e-/m^3)
% oxden = O+ density (O+/m^3)
% hyden = H+ density (H+/m^3)
% heden = He+ density (He+/m^3)
% hytot = H+ total tube content (H+)
% hetot = He+ total tube content (H+)
% lshv = L-shell vector of results

%Check all input are column vectors
mlt=mlt(:);
timmat=timmat(:);
dno=dno(:);
year=year(:);

re=6371;                                %Earth radi

%% Load in plasimod data
[lshcp,lshmx,lshv,mltv]=plasires2(dno(1),year(1));

%% Use interpolation to determine the H+ tube content and H+ and O+
%% densities along the required radial ray paths for points above
%% 350km.
hytot=zeros(length(timmat),length(lshv));
hetot=hytot;
hyden=hytot;
heden=hytot;
oxden=hytot;

for i=1:length(timmat)
 [hytotj,hetotj,hyeqj,heeqj,oxeqj,iricf,hes]=plasires3(timmat(i),dno(i),year(i),mlt(i)*ones(1,length(lshv)),lshv(:)');  
 hytot(i,:)=hytotj(:)';
 hyden(i,:)=hyeqj(:)';
 oxden(i,:)=oxeqj(:)';
 if hes==1
  hetot(i,:)=hetotj(:)';
  heden(i,:)=heeqj(:)';
 else
  hetot(i,:)=nan;
  heden(i,:)=nan;
 end
end 

%Electron density
if hes==1
 eden=hyden+heden+oxden;
else
 eden=hyden+oxden;
end

