%% msiscal.m calculates the neutral temperature, neutral oxygen, hydrogen,
%% helium and molecular nitrogen densities at the height of interest
%% using the inputted constants from the MSIS model.
%% Equations are from MSIS model (Hedin 1987 JGR)
%% 
%% Note - if hes==0 then hent and n2nt are not calculated and empty arrays
%%        are returned for these variables instead.
%%

function [tmpn,oxnt,hynt,hent,n2nt]=msiscal(tmpnex,tmpnb,tmpnsp,oxntb,hyntb,hentb,n2ntb,disb,dismat,hes)

% tmpnex = exospheric neutral temperature (K)
% tmpnb = base neutral temperature (K)
% tmpnsp = shaping constant
% oxntb = base neutral oxygen density (m^-3)
% hyntb = base neutral hydrogen density (m^-3)
% hentb = base neutral helium density (m^-3)
% n2ntb = base neutral molecular nitrogen density (m^-3)
% disb = base radial distance (km)
% dismat = radial distance of point of interest (km)
% hes = He+ simulated (0 = no, 1 = yes) (scalar)

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

ep=(dismat-disb).*disb./dismat;                                  %eqn A4a
tmpn=tmpnex-(tmpnex-tmpnb).*exp(-tmpnsp.*ep);                    %eqn A1a - neutral temperature
g0=9.8./(disb./re).^2;                                           %eqn A17a - base 350km gravity (m/s^2)
shhy=1e-3*(k.*tmpnex)./(amu*g0);                                 %neutral hydrogen scale height (km)
gm2hy=1./(shhy.*tmpnsp);                                         %eqn A16a - hydrogen gamma_2 
hynt=hyntb.*(tmpnb./tmpn).^(0.6+gm2hy).*exp(-ep./shhy);          %eqn A13, A14a & A15 - neutral hydrogen density
oxnt=oxntb.*(tmpnb./tmpn).^(1+16*gm2hy).*exp(-16*ep./shhy);      %eqn A13, A14a & A15 - neutral oxygen density
if hes==1
 hent=hentb.*(tmpnb./tmpn).^(0.6+4*gm2hy).*exp(-4*ep./shhy);     %eqn A13, A14a & A15 - neutral helium density
 n2nt=n2ntb.*(tmpnb./tmpn).^(1+28*gm2hy).*exp(-28*ep./shhy);     %eqn A13, A14a & A15 - neutral N2 density
else
 hent=[];
 n2nt=[];
end
