%% plasiray12.m is a modified version of plasipt15.m which allows allows
%% different time and date information for each field line. It calculates
%% the ion densities for points above 350km and outside the plasmasphere.
%%

function [oxdenp,hydenp,hedenp]=plasiray12(hyeq,heeq,lsh,latmatmp,dismatp,lonlsh,timmat,year,month,day,lshcp,nop,iricf,gms,hes,sms)

% hyeq = equatorial or end of field line (at distance lshcp) hydrogen ion density (H+/m^-3)
% heeq = equatorial or end of field line (at distance lshcp) helium ion density (He+/m^-3)
% lsh = L-shells of points of interest
% latmatmp = magnetic latitude of point of interest (rads)
% dismatp = radial distance of point of interest (km)
% lonlsh = magnetic longitude of point of interest (rads) 
% timmat = univerisal time of point of interest (fractional hours)
% year = year
% month = month
% day = day 
% lshcp = upper radial L-shell (scalar)
% nop = number of points per field line (scalar)
% iricf = IRI2001 FoF2 data set to be used (0 = CCIR, other = URSI) (scalar)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (scalar)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% oxdenp = O+ density at point of interest (m^-3)
% hydenp = H+ density at point of interest (m^-3)

%Make sure all input are column vectors
hyeq=hyeq(:);
heeq=heeq(:);
lsh=lsh(:);
latmatmp=latmatmp(:);
dismatp=dismatp(:);
lonlsh=lonlsh(:);
timmat=timmat(:);
year=year(:);
month=month(:);
day=day(:);

re=6371;                 %Earth's radius
r0=re+350;               %Base radial distance
rb=re+700;               %upper radial distance for IRI2001 to be used
idhf2=50;                %height above HmF2 at which use diffusive equilibrium for O+ (km)
rats=0.1;                %diffusive/chemical loss ratio at which to start upwards diffusive profile
sz=length(hyeq);         %number of points of interest

%Calculate the required solar parameters
dno=daynum(day,month,year);             %day number
[kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timmat,dno,year);

%Calculate the co-ordinate parameters of the field line
svy=ones(size(lsh));
latmb=lshpar(lsh,-svy,r0*svy);
latmb(latmatmp<0)=-latmb(latmatmp<0);
[lonmat,latmat,latmatm,dismat,dflm]=flpar3(latmb,lonlsh,r0,lshcp,nop);
sv=ones(1,nop);
svm=ones(size(dismat));

%Calculate required ionosphere electron density values at the base of the field line.
%This is done for every point between r0 and 700km altitude.
ionden=zeros(size(lonmat));             %electron density
hmf2s=ionden;                           %dummy matrix to hold HmF2 data
poi=dismat<rb;                          %points for which IRI2001 to be used
[pyi,pxi]=find(poi);                    %y and x index of points where IRI2001 to be used
svp=ones(sum(sum(poi)),1);
[ionden(poi),nmf2,hmf2s(poi)]=modiri5(lonmat(poi),latmat(poi),timmat(pyi),day(pyi),month(pyi),year(pyi),dismat(poi),rb,iricf,gms);
hmf2p=hmf2s(:,1);                       %HmF2 at base of field line
clear hmf2s nmf2

%Calculate the required neutral densities. 
apmat2=zeros(size(apmat,1),nop,7);
for i=1:nop
 apmat2(:,i,:)=apmat;
end
[oxnt,hynt,hent,n2nt,tmpn]=modmsis1(lonmat,latmat,dismat,timmat*sv,day*sv,month*sv,year*sv,f107a*sv,f107p*sv,apmat2,sms);

%Calculate the temperatures at the points of interest.
[tmpi,tmpe,tmpic,tmpec,tdiur,tdiurc]=temppaw2(lonmat,latmat,dismat,timmat*sv,year*sv,month*sv,day*sv,tmpn,tmpn,kp*sv,f107*sv);
tmpt=tmpi+tmpe;                                             %total temperature

%Determine the gravity component along the field at each point of interest.
g0=gravmod(dismat,latmatm);

%Set the oxygen density equal to the IRI2001 electron density values. It is assumed
%that the oxygen density is equal to the electron density for all points below
%HmF2+idhf2 km (using the HmF2 value at the base of the field line). For points 
%above this diffusive equilibrium will be used.
oxden=ionden;                           %oxygen density at points of interest
shoxm=0*svm;                            %O+ scale heights

%Take initial H+ and He+ densities equal to the polar densities
[hydeni,hedeni]=polarmod(latmatm,dismat,kp);
if hes==0
 hedeni=0*hedeni;
end

%Use the diffusive equilibrium model to calculate the O+ density profiles.
[oxden,shoxm]=plasipt13(oxden,0,shoxm,hydeni,hedeni,lsh,dismat,dflm,latmatm,g0,tmpi,tmpe,tmpt,(re+hmf2p+idhf2),1);

%Use the diffusive equilibrium model to calculate the He+ density profiles. 
if hes==1
 heden=plasipt16(oxden,hydeni,heeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp);
else
 heden=hedeni;
end

%Use the diffusive equilibrium model to calculate the H+ density profiles. 
hyden=plasipt14(oxden,heden,hyeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp);

%Recalculate the O+ density profiles with H+ and He+ density profiles.
oxden=plasipt13(oxden,0,shoxm,hyden,heden,lsh,dismat,dflm,latmatm,g0,tmpi,tmpe,tmpt,(re+hmf2p+idhf2),0);

%Recalculate the He+ density profiles with O+ and H+ density profiles. 
if hes==1
 [heden,shhem]=plasipt16(oxden,hyden,heeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp);
else
 shhem=heden;
end

%Recalculate the H+ density profiles with O+ density profiles. 
[hyden,shhym]=plasipt14(oxden,heden,hyeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp);

%Modifiy the H+ and He+ densities if required to include the effects of chemical equilibrium. 
hyden=plsd4(hyden,oxden,shhym,shoxm,dismat,dflm,tmpi,tmpe,tmpt,tmpn,oxnt,hynt,g0,r0,rats);
if hes==1
 heden=plsd6(heden,oxden,shhem,shoxm,dismat,dflm,tmpi,tmpe,tmpt,tmpn,n2nt,hent,g0,tdiur,f107,f107a,r0,rats);
end

%Use interpolation to determine the ion densities
dflmp=distfl2(latmatmp,lsh,r0*ones(sz,1));   %distance point of interest along field line
oxdenp=intfl(dflm,oxden,dflmp,0);            %O+ density at point of interest
hydenp=intfl(dflm,hyden,dflmp,0);            %H+ density at point of interest
if hes==1
 hedenp=intfl(dflm,heden,dflmp,0);            %He+ density at point of interest
else
 hedenp=nan*ones(size(hydenp));
end
 

