%% plasieq1.m determines the electron densities for the specified magnetic
%% equatorial 'slice' at the given UT, day, month and year using data files
%% produced by plasimod(2).m. The results returned are from the closest UT
%% requested.
%%
%% NOTE - while the electron densities are returned are returned for all
%%        points, the oxygen and hydrogen densities are only returned above
%%        350km (zero is returned below 350km). 
%%
%% NOTE - minimum allowed altitude is 1000km, maximum altitude set by lshcp
%%        value from plasimod(2) data sets.
%%
%% NOTE - altitude term approximate, the closest actual values above maximum
%%        from plasimod(2) returned.
%%
%% 10/11/2001 - updated 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 [eden,oxden,hyden,heden,lonmatm,dismat,mlt,hes]=plasieq1(altmax,timmat,dno,year)

% altmax = maximum altitude of each radial path (km) (scalar)
% timmat = univerisal time (fractional hours) (scalar)
% dno = day number (scalar)
% year = year (scalar)
%
% eden = electron density (e-/m^3)
% oxden = O+ density (O+/m^3)
% hyden = H+ density (H+/m^3)
% heden = He+ density (He+/m^3)
% lonmatm = magnetic longitude of grid
% dismat = radial distance of grid
% mlt = magnetic local time vector
% hes = He+ simulated (0 = no, 1 = yes) (scalar)

re=6371;                                %Earth radi

%% Load in plasimod data set for required day
[hytotil,hetotil,hyeqil,heeqil,oxeqil,hytotiu,hetotiu,hyeqiu,heeqiu,oxeqiu,timmatl,timmatu,lsh,mlt,lshcp,gms,iricf,hes]=plasires(timmat,dno,year);

tim1=timmat-timmatl;
tim2=timmatu-timmat;
if tim1<tim2
 hytoti=hytotil;
 hetoti=hetotil;
 hyeqi=hyeqil;
 heeqi=heeqil;
 oxeqi=oxeqil;
else 
 hytoti=hytotiu;
 hetoti=hetotiu;
 hyeqi=hyeqiu;
 heeqi=heeqiu;
 oxeqi=oxeqiu;
end 
 
%% Select out the required components of data files that correspond to given
%% altitude minimums and maximums.
altv=re*(lsh-1);                   %equatorial altitude vector
poi=sum(altv<altmax)+1;            %index of altitude band above maximum altitude
hyden1=hyeqi(2:poi,:);             %required H+ densities from plasimod(2)
if hes==1
 heden1=heeqi(2:poi,:);            %required He+ densities from plasimod(2)
end
oxden1=oxeqi(2:poi,:);             %required O+ densities from plasimod(2)
if hes==1
 eden1=hyden1+heden1+oxden1;       %electron densities
else
 eden1=hyden1+oxden1;              %electron densities
end

%% Generate electron densities only for region between 90km and 350km.
disv=re+(90:26:350)';                                            %displacement vector
svx=ones(1,length(mlt));
svy=ones(length(disv),1);
svm=ones(length(disv),length(mlt));
lonmatmv=timml2lon(timmat*svx,mlt,dno*svx);                      %magnetic longitude vector
[lonmatv,latmatv]=mag2geo3(lonmatmv,0*svx);                      %geographic co-ordinates 
[day,month]=daynum2(dno,year);                                   %day and month
eden2=modiri5(svy*lonmatv,svy*latmatv,timmat*svm,day*svm,month*svm,year*svm,disv*svx,500,iricf,gms);

%% Produce the lookup tables required for output
eden=[eden2;eden1];                     %electron densities (e-/m^3)
hyden2=zeros(size(eden2));               
hyden=[hyden2;hyden1];                  %H+ densities (H+/m^3)
if hes==1
 heden=[hyden2;heden1];                 %He+ densities (H+/m^3)
end
oxden=[hyden2;oxden1];                  %O+ densities (O+/m^3)
lonmatm=ones(size(eden,1),1)*lonmatmv;  %magnetic longitude (rads)
dismat=[disv*svx;(re*lsh(2:poi))*svx];  %radial displacement (km)

if hes==0
 heden=nan*ones(size(hyden));
end