%% msisap8.m estimates the neutral oxygen and hydrogen 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 msisap8.m should come from msisap9.m, since it requires
%% the same L-shells to be used at each magnetic longitude.
%%
%% (21/6/00) It has be found that the simple approximation of assuming
%% vertical field lines along which to calculate the neutral parameters
%% can result in derivations of up to several hundred percent at 
%% altitudes around ~1500km. Therefore replace this approach by using
%% interpolation to determine the base values in each parameter at
%% 350km (using claculated values from msisap9.m) at the magnetic
%% latitude of interest and then use the MSIS equations to determine 
%% the neutral temperature and densities at each point along each 
%% field line.
%% 

function [oxntm,hyntm,tmpnm,oxntmc,hyntmc,tmpnmc]=msisap8(lshv,lshcp,timi,dismat,latmatm,oxnt350,hynt350,tmpn350,tmpnex,tmpnsp,oxnt350c,hynt350c,tmpn350c,tmpnexc,tmpnspc)

% 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)
% oxnt350c = conjugate neutral oxygen density at 350 km (m^-3)
% hynt350 = neutral hydrogen density at 350 km (m^-3)
% hynt350c = conjuagte neutral hydrogen density at 350 km (m^-3)
% tmpn350 = neutral temperature at 350 km (K)
% tmpn350c = conjugate neutral temperature at 350 km (K)
% tmpnex = neutral exospheric temperature (K)
% tmpnexc = conjuagte neutral exospheric temperature (K)
% tmpnsp = neutral temperature shaping parameter
% tmpnspc = conjugate neutral temperature shaping parameter
%
% oxntm = neutral oxygen density matrix (m^-3)
% hyntm = neutral hydrogen 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)
% 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);                   %required oxygen density at 350km 
hy35=hynt350(:,timi);                   %hydrogen density at 350km
tn35=tmpn350(:,timi);                   %neutral temperature at 350km
tnex=tmpnex(:,timi);                    %neutral exospheric temperature 
tnsp=tmpnsp(:,timi);                    %neutral temperature shaping parameter
ox35c=oxnt350c(:,timi);                 %required conjugate oxygen density at 350km 
hy35c=hynt350c(:,timi);                 %hydrogen conjugate density at 350km
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;
tmpnm=oxntm;
oxntmc=oxntm;
hyntmc=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
 
 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)];
 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));

 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));

 [tmpni,oxnti,hynti]=msiscal(tnexi,tnbi,tnspi,ox35i,hy35i,dismat(poi,1)*ones(1,sz(2)),dismat(poi,:));
 oxnti(oxnti<1)=1;
 oxntm(poi,:)=oxnti;
 hyntm(poi,:)=hynti;
 tmpnm(poi,:)=tmpni;

 [tmpnic,oxntic,hyntic]=msiscal(tnexic,tnbic,tnspic,ox35ic,hy35ic,dismat(poi,1)*ones(1,sz(2)),dismat(poi,:));
 oxntic(oxntic<1)=1;
 oxntmc(poi,:)=oxntic;
 hyntmc(poi,:)=hyntic;
 tmpnmc(poi,:)=tmpnic;

end
