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

function [oxdenp,hydenp,hedenp]=plasiray11(hyeq,heeq,oxeq,lsh,lonm,latm,timmat,day,month,year,nop,iricf,gms,hes,sms)

% hyeq = equatorial H+ density (H+/m^3)
% heeq = equatorial He+ density (He+/m^3)
% oxeq = equatorial O+ density (H+/m^3)
% lsh = L-shell of points of interest
% lonm = magnetic longitude of points of interest
% latm = magnetic latitude of points of interest
% timmat = univerisal time (fractional hours)
% day = day of points of interest
% month = month of points of interest
% year = year of points of interest
% 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)
%
% oxenp = O+ density of point of interest (O+/m^3)
% hyenp = H+ density of point of interest (H+/m^3)
% heenp = He+ density of point of interest (He+/m^3)

%Make sure all input are columns
hyeq=hyeq(:);
heeq=heeq(:);
oxeq=oxeq(:);
lsh=lsh(:);
lonm=lonm(:);
latm=latm(:);
timmat=timmat(:);
day=day(:);
month=month(:);
year=year(:);

re=6371;                 %Earth radi (km)
r0=re+350;               %base radial distance of field lines
rb=re+700;               %upper radial distance for IRI95 to be used
sz=length(hyeq);         %number of points of interest
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

%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
[lonmat2,lonmat2c,latmat2,latmat2c,latmatm,dismat,dflm]=flpar(lsh,lonm,r0,nop);
sv=ones(1,nop);
svm=ones(size(dismat));

%Because only half the field line is being considered, need to swap to conjugate
%geographic co-ordinates for points in southern magnetic hemisphere
lonmat=lonmat2;
latmat=latmat2;
pois=latm<0;
if any(pois)
 lonmat(pois,:)=lonmat2c(pois,:);
 latmat(pois,:)=latmat2c(pois,:);
end 

%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 IRI95 to be used
[pyi,pxi]=find(poi);                    %y and x index of points where IRI95 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=plasipt13(oxden,oxeq,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,shoxm]=plasipt13(oxden,oxeq,shoxm,hyden,heden,lsh,dismat,dflm,latmatm,g0,tmpi,tmpe,tmpt,(re+hmf2p+idhf2),1);

%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(latm,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


