%% modmsis4.m is an improved version of msisap8.m that estimates the neutral 
%% oxygen, hydrogen, helium and molecular N2 densities, and the neutral 
%% temperature at the points of interest and the conjugate points. The 
%% inputted values correspond to the required parameters at 350km, from which
%% the equations from MSIS-90 are used to calculate the required parameters. 
%% The equations used from MSIS-90 are listed in the code and are from the 
%% 1987 Hedin paper in May JGR (Vol 92 A5 p4649). Note that they are modifed 
%% to use 350km as their base height, not 120km.
%% The input for modmsis4.m should come from modmsis3.m, since it requires
%% the same L-shells to be used at each magnetic longitude.
%%
%% 6/3/2002 - include option to not calculate neutral parameters required by He+
%%            if He+ is not to be simulated.
%%

function [oxntm,hyntm,hentm,n2ntm,tmpnm,oxntmc,hyntmc,hentmc,n2ntmc,tmpnmc]=modmsis4(lshv,lshcp,timi,dismat,latmatm,oxnt350,hynt350,hent350,n2nt350,tmpn350,tmpnex,tmpnsp,oxnt350c,hynt350c,hent350c,n2nt350c,tmpn350c,tmpnexc,tmpnspc,hes)

% lshv = L-shell vector used at each magnetic longitude
% lshcp = maximum L-shell for which field line calculations are undertaken (scalar)
% timi = time index to be considered (scalar)
% dismat = radial distance of points of interests along field lines
% latmatm = magnetic latitude of points of interest along field line
% oxnt350 = neutral oxygen density at 350 km (m^-3)
% hynt350 = neutral hydrogen density at 350 km (m^-3)
% hent350 = neutral helium density at 350 km (m^-3)
% n2nt350 = neutral N2 density at 350 km (m^-3)
% tmpn350 = neutral temperature at 350 km (K)
% tmpnex = neutral exospheric temperature (K)
% tmpnsp = neutral temperature shaping parameter
% oxnt350c = conjugate neutral oxygen density at 350 km (m^-3)
% hynt350c = conjugate neutral hydrogen density at 350 km (m^-3)
% hent350c = conjugate neutral helium density at 350 km (m^-3)
% n2nt350c = conjugate neutral N2 density at 350 km (m^-3)
% tmpn350c = conjugate neutral temperature at 350 km (K)
% tmpnexc = conjugate neutral exospheric temperature (K)
% tmpnspc = conjugate neutral temperature shaping parameter
% hes = He+ simulated (0 = no, 1 = yes) (scalar)
%
% oxntm = neutral oxygen density matrix (m^-3)
% hyntm = neutral hydrogen density matrix (m^-3)
% hentm = neutral helium density matrix (m^-3)
% n2ntm = neutral N2 density matrix (m^-3)
% tmpnm = neutral temperature matrix (K)
% oxntmc = conjugate neutral oxygen density matrix (m^-3)
% hyntmc = conjugate neutral hydrogen density matrix (m^-3)
% hentmc = conjugate neutral helium density matrix (m^-3)
% n2ntmc = conjugate neutral N2 density matrix (m^-3)
% tmpnmc = conjugate neutral temperature matrix (K)

re=6371;                                %Earth radi (km)
k=1.38e-23;                             %Boltzman constant (J/K) 
amu=1.66e-27;                           %atomic mass unit (kg)

ox35=oxnt350(:,timi);                   %oxygen density at 350km 
hy35=hynt350(:,timi);                   %hydrogen density at 350km
if hes==1
 he35=hent350(:,timi);                  %helium density at 350km
 n235=n2nt350(:,timi);                  %N2 density at 350km
end
tn35=tmpn350(:,timi);                   %neutral temperature at 350km
tnex=tmpnex(:,timi);                    %neutral exospheric temperature 
tnsp=tmpnsp(:,timi);                    %neutral temperature shaping parameter
ox35c=oxnt350c(:,timi);                 %conjugate oxygen density at 350km 
hy35c=hynt350c(:,timi);                 %conjugate hydrogen density at 350km
if hes==1
 he35c=hent350c(:,timi);                %conjugate helium density at 350km
 n235c=n2nt350c(:,timi);                %conjugate N2 density at 350km
end
tn35c=tmpn350c(:,timi);                 %neutral conjugate temperature at 350km
tnexc=tmpnexc(:,timi);                  %neutral conjugate exospheric temperature 
tnspc=tmpnspc(:,timi);                  %neutral conjugate temperature shaping parameter

%Determine required L-shells for each magnetic longitude
lsh=lshv;
lsh(lsh>lshcp)=[];

sz=size(dismat);
nopy=length(lsh);                       %number of field lines per magnetic longitude
nopx=sz(1)/nopy;                        %number of magnetic longitudes

%Result matrices
oxntm=zeros(sz);
hyntm=oxntm;
hentm=oxntm;
n2ntm=oxntm;
tmpnm=oxntm;
oxntmc=oxntm;
hyntmc=oxntm;
hentmc=oxntm;
n2ntmc=oxntm;
tmpnmc=oxntm;

%Loop through each magnetic longitude to determine the required parameters
for i=1:nopx
 
 id1=1+(i-1)*nopy;                 %lower i'th index 
 id2=i*nopy;                       %upper i'th index 
 poi=(id1:id2)';                   %i'th points of interest
 
 if hes==1
  dat1=[flipud(tnexc(poi,1)) flipud(tn35c(poi,1)) flipud(tnspc(poi,1)) flipud(ox35c(poi,1)) flipud(hy35c(poi,1)) flipud(he35c(poi,1)) flipud(n235c(poi,1))];
  dat2=[tnex(poi,1) tn35(poi,1) tnsp(poi,1) ox35(poi,1) hy35(poi,1) he35(poi,1) n235(poi,1)];
 else
  dat1=[flipud(tnexc(poi,1)) flipud(tn35c(poi,1)) flipud(tnspc(poi,1)) flipud(ox35c(poi,1)) flipud(hy35c(poi,1))];
  dat2=[tnex(poi,1) tn35(poi,1) tnsp(poi,1) ox35(poi,1) hy35(poi,1)];
 end
 latmi=latmatm(poi,:);

 resm=interp1([-flipud(latmatm(poi,1));latmatm(poi,1)],[dat1;dat2],latmi(:));
  
 tnexi=reshape(resm(:,1),nopy,sz(2));
 tnbi=reshape(resm(:,2),nopy,sz(2));
 tnspi=reshape(resm(:,3),nopy,sz(2));
 ox35i=reshape(resm(:,4),nopy,sz(2));
 hy35i=reshape(resm(:,5),nopy,sz(2));
 if hes==1
  he35i=reshape(resm(:,6),nopy,sz(2));
  n235i=reshape(resm(:,7),nopy,sz(2));
 else
  he35i=[];
  n235i=[];
 end

 resmc=interp1([-flipud(latmatm(poi,1));latmatm(poi,1)],[dat1;dat2],-latmi(:));

 tnexic=reshape(resmc(:,1),nopy,sz(2));
 tnbic=reshape(resmc(:,2),nopy,sz(2));
 tnspic=reshape(resmc(:,3),nopy,sz(2));
 ox35ic=reshape(resmc(:,4),nopy,sz(2));
 hy35ic=reshape(resmc(:,5),nopy,sz(2));
 if hes==1
  he35ic=reshape(resmc(:,6),nopy,sz(2));
  n235ic=reshape(resmc(:,7),nopy,sz(2));
 else
  he35ic=[];
  n235ic=[];
 end

 [tmpni,oxnti,hynti,henti,n2nti]=msiscal(tnexi,tnbi,tnspi,ox35i,hy35i,he35i,n235i,dismat(poi,1)*ones(1,sz(2)),dismat(poi,:),hes);

 oxnti(oxnti<1)=1;
 n2nti(n2nti<1)=1;
 oxntm(poi,:)=oxnti;
 hyntm(poi,:)=hynti;
 if hes==1
  hentm(poi,:)=henti;
  n2ntm(poi,:)=n2nti;
 end
 tmpnm(poi,:)=tmpni;

 [tmpnic,oxntic,hyntic,hentic,n2ntic]=msiscal(tnexic,tnbic,tnspic,ox35ic,hy35ic,he35ic,n235ic,dismat(poi,1)*ones(1,sz(2)),dismat(poi,:),hes);

 oxntic(oxntic<1)=1;
 n2ntic(n2ntic<1)=1;
 oxntmc(poi,:)=oxntic;
 hyntmc(poi,:)=hyntic;
 if hes==1
  hentmc(poi,:)=hentic;
  n2ntmc(poi,:)=n2ntic;
 end
 tmpnmc(poi,:)=tmpnic;
end

if hes==0
 hentm=[];
 n2ntm=[];
 hentmc=[];
 n2ntmc=[];
end