%% plasiray1.m calculates the electron, and ion densities along the
%% ray paths specified by the user using the data tables produced
%% by plasimod(2). For each ray path it also determines the TEC, 
%% mean ionospheric height. All inputs should be vectors of the same
%% size, while the outputs are of the format where each row corresponds
%% to a ray path, except for tec, mih, lonmih and latmih which are column
%% vectors. Note that mih (mean ionospheric height) is really a misnomer
%% since the ray specificed by the user does not have to pass through
%% the ionosphere at all. It actually represents the point along the given
%% ray path where half the TEC is below and half is above.
%% NOTE: Points between 0 and 100km have 0 returned for all densities.
%%
%% 14/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 [tec,mih,lonmih,latmih,eden,oxden,hyden,heden,lonmat,latmat,dismat]=plasiray1(loni,lonf,lati,latf,alti,altf,noprp,lt,timmat,day,month,year,nop)

% loni = initial geographic longitude vector of ray paths (rads)
% lonf = final geographic longitude vector of ray paths (rads)
% lati = initial geographic latitude vector of ray paths (rads)
% latf = final geographic latitude vector of ray paths (rads)
% alti = initial altitude of ray paths (km)
% altf = final altitude of ray paths (km)
% noprp = number of points along ray paths (scalar)
% lt = line spacing interval type (0 = log, other for linear) (scalar)
% timmat = UT vector
% day = day vector
% month = month vector
% year = year vector
% nop = number of points per field line (scalar)
%
% tec = TEC vector 
% mih = 'mean ionospheric height' vector (km)
% lonmih = 'mean ionospheric height' longitude vector (rads)
% latmih = 'mean ionospheric height' latitude vector (rads)
% eden = electron densities (e-/m^3)
% oxden = O+ densities (O+/m^3)
% hyden = H+ densities (H+/m^3)
% heden = He+ densities (He+/m^3)
% lonmat = geographic longitude (rads)
% latmat = geographic latitude (rads)
% dismat = radial distance (km)


%Make sure all input are column vectors
loni=loni(:);
lonf=lonf(:);
lati=lati(:);
latf=latf(:);
alti=alti(:);
altf=altf(:);
timmat=timmat(:);
day=day(:);
month=month(:);
year=year(:);

re=6371;                      %Earth radi (km)

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

%Determine lshcp parameter
lshcp=plasires2(dno,year);

%% Convert from spherical to carteian co-ordinates to determine the true distance
%% between the start and end of each eay path
[xi,yi,zi]=sph2cart(loni,lati,(alti+re));         %initial points cartesian co-ords
[xf,yf,zf]=sph2cart(lonf,latf,(altf+re));         %final points cartesian co-ords
dist=sqrt((xf-xi).^2+(yf-yi).^2+(zf-zi).^2);      %total length of ray paths (km)
distx=xf-xi;                                      %x component distance (km)
disty=yf-yi;                                      %y component distance (km)
distz=zf-zi;                                      %z component distance (km)

if lt==0                                          %log scaling
 dv=logspace(-2,0,(noprp-1));                     %fractional vector to break field line up into
 dsf=[dv-dv(1) 1];                                %log distance scaling factor
else                                              %linear scaling 
 dsf=linspace(0,1,noprp);                         %linear distance scaling factor
end

sv=ones(1,noprp);     
xco=xi*sv+distx*dsf;                              %required x-coordinates
yco=yi*sv+disty*dsf;                              %required y-coordinates
zco=zi*sv+distz*dsf;                              %required z-coordinates

%% Convert back from cartesian to spherical co-ordinates
[lonmat,latmat,dismat]=cart2sph(xco,yco,zco);

%% Reset 1st and last radial displacements to agree with input (stops round off error)
dismat(:,1)=re+alti;
dismat(:,noprp)=re+altf;

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

%% For points with altitudes less than 0 km set densities to NaN. 
poie=dismat<re;
if any(any(poie))
 eden(poie)=nan;
 oxden(poie)=nan;
 hyden(poie)=nan;
 heden(poie)=nan;
end 

%% Calculate the electron densities between 100km and 350km using the IRI2001
poi=dismat>=(re+100) & dismat<=(re+350);          %points between 100km and 350km
if any(any(poi))
 [pyi,pxi]=find(poi);
 eden(poi)=modiri5(lonmat(poi),latmat(poi),timmat(pyi),day(pyi),month(pyi),year(pyi),dismat(poi),500,iricf,gms);
end

%% For points between 0km and 100km altitude, assume zero density.
poi100=dismat>=re & dismat<(re+100);

%% Calculate the L-shell of each point of interest above 350km. Points with L-shell
%% less than min(7,lshcp) are considered to be in the plasmasphere while points 
%% outside are 'polar'. 
svm=ones(size(dismat));
[lonmatm,latmatm]=geo2mag3(lonmat,latmat);             %magnetic co-ordinates
lshm=lshpar(-svm,latmatm,dismat);                      %L-shell of points of interest
poips=~poie & ~ poi & ~poi100 & lshm<=min(lshcp,7);    %'plasmasphere' points
poipo=~poie & ~ poi & ~poi100 & ~poips;                %'polar' points

%% 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.
hyeqm=0*svm;                                           %H+ equatorial results matrix
heeqm=hyeqm;                                           %He+ equatorial results matrix
oxeqm=hyeqm;                                           %O+ equatorial results matrix
svx=ones(1,size(lonmatm,2));
mltm=timut2ml2(timmat*svx,lonmatm,dno*svx);            %MLT matrix

for j=1:length(timmat);                                %loop through each ray path
 [hytot,hetot,hyeqp,heeqp,oxeqp]=plasires3(timmat(j),dno(j),year(j),mltm(j,:),lshm(j,:));
 hyeqm(j,:)=hyeqp;                                     %record j'th H+ equatorial densities
 oxeqm(j,:)=oxeqp;                                     %record j'th O+ equatorial densities
 if hes==1
  heeqm(j,:)=heeqp;                                    %record j'th He+ equatorial densities
 else
  heeqm(j,:)=nan;
 end
end

%% Calculate the required densities for points inside the plasmasphere.
if any(any(poips))
 [pyi,pxi]=find(poips);
 [oxden(poips),hyden(poips),heden(poips)]=plasiray11(hyeqm(poips),heeqm(poips),oxeqm(poips),lshm(poips),lonmatm(poips),latmatm(poips),timmat(pyi),day(pyi),month(pyi),year(pyi),nop,iricf,gms,hes,sms);
 if hes==1
  eden(poips)=oxden(poips)+hyden(poips)+heden(poips);
 else
  eden(poips)=oxden(poips)+hyden(poips);
 end
end

%% Calculate the required densities for points in the polar regions.
if any(any(poipo))
 [pyi,pxi]=find(poipo);
 [oxden(poipo),hyden(poipo),heden(poipo)]=plasiray12(hyeqm(poipo),heeqm(poipo),lshm(poipo),latmatm(poipo),dismat(poipo),lonmatm(poipo),timmat(pyi),year(pyi),month(pyi),day(pyi),lshcp,nop,iricf,gms,hes,sms);
 if hes==1
  eden(poipo)=oxden(poipo)+hyden(poipo)+heden(poipo);
 else
  eden(poipo)=oxden(poipo)+hyden(poipo);
 end
end 

%% Calculate the required mean ionospheric height parameters
tec=zeros(length(loni),1);              %TEC results vector
mih=tec;                                %mean ionospheric height vector
lonmih=tec;                             %geogrpahic longitude of mean ionospheric height
latmih=tec;                             %geogrpahic latitude of mean ionospheric height
pnan=sum(isnan(eden),2);                %check for NaN
poi=pnan==0;                            %terms to calculate ionospheric parameters
if any(poi)
 [tec(poi),mih(poi),lonmih(poi),latmih(poi)]=teccal3(eden(poi,:),dismat(poi,:),lonmat(poi,:),latmat(poi,:),1);
end

