%% plasisc12.m is designed to be called by plasisc11.m. It calculates the
%% oxygen, hydrogen and helium ion densities profiles at selected heights 
%% along the given L-shell. 
%%
%% 16/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 [hyden,heden,oxden,hydenc,hedenc,oxdenc,latmatm,dismat]=plasisc12(hyeq,heeq,oxeq,lsh,lonlsh,timmat,year,month,day,nop,iricf,gms,hes,sms)

% hyeq = equatorial hydrogen ion density of flux tube (m^-3)
% heeq = equatorial helium ion density of flux tube (m^-3)
% oxeq = equatorial oxygen ion density of flux tube (m^-3)
% lsh = L-shells of interest
% lonlsh = magnetic longitude of L-shells (rads) (scalar or matrix size(hyeq))
% timmat = univerisal time (fractional hours) (scalar)
% year = year (scalar)
% month = month (scalar)
% day = day (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)
%
% hyden = northern hemisphere hydrogen density profile (m^-3)
% heden = northern hemisphere helium density profile (m^-3)
% oxden = northern hemisphere oxygen density profile (m^-3)
% hydenc = southern hemisphere hydrogen density profile (m^-3)
% hedenc = southern hemisphere helium density profile (m^-3)
% oxdenc = southern hemisphere oxygen density profile (m^-3)
% latmatm = magnetic latitude of profile points (all positive)
% dismat = radial displacement of profile points (km)

re=6371;                                %Earth's radius
r0=re+350;                              %Base radial distance
rb=re+700;                              %upper radial distance for IRI95 to be used
rats=0.1;                               %diffusive/chemcial ratio from which to start upwards diffusive profile
idhf2=50;                               %height above HmF2 at which use diffusive equilibrium for O+ (km)

%convert all input to column vectors
hyeq=hyeq(:);
heeq=heeq(:);
oxeq=oxeq(:);
lsh=lsh(:);
lonlsh=lonlsh(:);

%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(hyeq));
if prod(size(lonlsh))==1
 lonlsh=lonlsh*svy;
end 
[lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm]=flpar(lsh,lonlsh,r0,nop);
sv=ones(1,size(lonmat,2));
svm=ones(size(lonmat));

%Calculate required ionosphere electron density values for the two bases of the field line.
%This is done for every point between r0 and 700km altitude along each field line.
ionden=zeros(size(lonmat));             %electron density
iondenc=ionden;                         %conjugate electron density
hmf2s=ionden;                           %dummy matrix to hold HmF2 data
poi=dismat<rb;                          %points for which IRI95 to be used
[py,px]=find(poi);                      %index of poi
sv2=[1 1];

loni=lonmat(poi);
lonic=lonmatc(poi);
lati=latmat(poi);
latic=latmatc(poi);
dism=dismat(poi);
sv3=ones(size(loni(:)))*sv2;
[idn,nmf2,hmf2]=modiri5([loni(:) lonic(:)],[lati(:) latic(:)],timmat*sv3,day*sv3,month*sv3,year*sv3,dism(:)*sv2,rb,iricf,gms);
ionden(poi)=idn(:,1);                   %record electron densities
iondenc(poi)=idn(:,2);                  %record conjugate electron densities
hmf2s(poi)=hmf2(:,1);                   %record HmF2 data
hmf2p=hmf2s(:,1);                       %HmF2 at base of field line
hmf2s(poi)=hmf2(:,2);                   %reocrd conjugate HmF2 data
hmf2pc=hmf2s(:,1);                      %conjugate HmF2 at base of field line

%Calculate the required neutral densities. 
sz=size(lonmat);
[oxnm,hynm,henm,n2nm,tmpnt]=modmsis1([lonmat(:) lonmatc(:)],[latmat(:) latmatc(:)],dismat(:)*sv2,timmat,day,month,year,f107a,f107p,apmat,sms);

oxnt=reshape(oxnm(:,1),sz);             %record neutral oxygen densities
hynt=reshape(hynm(:,1),sz);             %record neutral hydrogen densities
hent=reshape(henm(:,1),sz);             %record neutral helium densities
n2nt=reshape(n2nm(:,1),sz);             %record neutral N2 densities
oxntc=reshape(oxnm(:,2),sz);            %record conjugate neutral oxygen densities
hyntc=reshape(hynm(:,2),sz);            %record conjugate neutral hydrogen densities
hentc=reshape(henm(:,2),sz);            %record conjugate neutral helium densities
n2ntc=reshape(n2nm(:,2),sz);            %record conjugate neutral N2 densities
tmpn=reshape(tmpnt(:,1),sz);            %record neutral temperature
tmpnc=reshape(tmpnt(:,2),sz);           %record conjugate neutral temperature

%Calculate the temperatures at the points of interest and conjugate points.
[tmpi,tmpe,tmpic,tmpec,tdiur,tdiurc]=temppaw2(lonmat,latmat,dismat,timmat*svm,year*svm,month*svm,day*svm,tmpn,tmpnc,kp*svm,f107*svm);
tmpt=tmpi+tmpe;                                             %total temperature
tmptc=tmpic+tmpec;                                          %total conjugate 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
oxdenc=iondenc;                         %conjugate 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);
oxdenc=plasipt13(oxdenc,oxeq,shoxm,hydeni,hedeni,lsh,dismat,dflm,latmatm,g0,tmpic,tmpec,tmptc,(re+hmf2pc+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);
 hedenc=plasipt16(oxdenc,hydeni,heeq,lsh,-latmatm,dismat,dflm,g0,tmpic,tmpec,tmptc,kp);
else
 heden=hedeni;
 hedenc=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);
hydenc=plasipt14(oxdenc,hedenc,hyeq,lsh,-latmatm,dismat,dflm,g0,tmpic,tmpec,tmptc,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);
[oxdenc,shoxmc]=plasipt13(oxdenc,oxeq,shoxm,hydenc,hedenc,lsh,dismat,dflm,latmatm,g0,tmpic,tmpec,tmptc,(re+hmf2pc+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);
 [hedenc,shhemc]=plasipt16(oxdenc,hydenc,heeq,lsh,-latmatm,dismat,dflm,g0,tmpic,tmpec,tmptc,kp);
else
 shhem=heden;
 shhemc=hedenc;
end

%Recalculate the H+ density profiles with O+ and He+ density profiles. 
[hyden,shhym]=plasipt14(oxden,heden,hyeq,lsh,latmatm,dismat,dflm,g0,tmpi,tmpe,tmpt,kp);
[hydenc,shhymc]=plasipt14(oxdenc,hedenc,hyeq,lsh,-latmatm,dismat,dflm,g0,tmpic,tmpec,tmptc,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);
hydenc=plsd4(hydenc,oxdenc,shhymc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,oxntc,hyntc,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);
 hedenc=plsd6(hedenc,oxdenc,shhemc,shoxmc,dismat,dflm,tmpic,tmpec,tmptc,tmpnc,n2ntc,hentc,g0,tdiurc,f107,f107a,r0,rats);
end
