%% plasipt1.m determines the electron densities at each of the specified
%% locations (longitude,latitude and altitude) all at the same UT, day,
%% month and year using data files produced by plasimod1.m. It is a
%% general purpose program within the above limitations, if the user
%% requires a magnetic longitude slice or a global magnetic electron
%% density calculation then plasisc.m and plasigb.m should be used instead
%% since they are specifically written to calculate this.
%%
%% NOTE - while the electron densities are returned for all
%%        points, the oxygen and hydrogen densities are only returned above
%%        350km (zero is returned below 350km). 
%%        
%% NOTE - valid altitudes are between 90km to a maximum given by the rounded
%%        down lshcp value from the plasmod(2) results. Values outside this
%%        range will return -1 for all output.
%% 
%% 12/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 [eden,oxden,hyden,heden]=plasipt1(lonmat,latmat,altmat,timmat,day,month,year,nop)

% lonmat = geographic longitude of points of interest (rads)
% latmat = geographic latitude of points of interest (rads)
% altmat = altitude of points of interest (km)
% timmat = univerisal time (fractional hours) (scalar)
% day = day (scalar)
% month = month (scalar)
% year = year (scalar)
% nop = number of points per field line (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)

re=6371;                                %Earth radi
dismat=altmat+re;                       %radial distance matrix

%Result matrices
eden=-ones(size(lonmat));               
oxden=eden;
hyden=eden;
heden=eden;

%Convert non-scalar input to column vectors
lonmat=lonmat(:);
latmat=latmat(:);
dismat=dismat(:);

%Determine parameters used in plasimod.m
dno=daynum(day,month,year);
[gms,iricf,hes,sms]=plasires4(dno,year);

%Determine the maximum allowed altitude.
dno=daynum(day,month,year);
lshcp=plasires2(dno,year);
altmx=1e3*floor(lshcp*re/1e3);

poi1=altmat>=90 & altmat<=350;          %points between 90km and 350km
poi2=altmat>350 & altmat<=altmx;        %points between 350km and maximum altitude

%% Calculate the electron densities between 90km and 350km using the IRI2001
if any(poi1)
 svi=ones(sum(poi1),1);
 eden(poi1)=modiri5(lonmat(poi1),latmat(poi1),timmat*svi,day*svi,month*svi,year*svi,dismat(poi1),500,iricf,gms);
end

%% Calculate the ion and electron densites between 350km and maximum altitude
if any(poi2)
 [oxden(poi2),hyden(poi2),heden(poi2)]=plasipt11(lonmat(poi2),latmat(poi2),dismat(poi2),timmat,day,month,year,lshcp,nop,iricf,gms,hes,sms);
 if hes==1
  eden(poi2)=oxden(poi2)+hyden(poi2)+heden(poi2);
 else
  eden(poi2)=oxden(poi2)+hyden(poi2);
 end
end
