%% modmsis3.m is an improved version of msisap9.m. It calls the MSIS model 
%% to calculate neutral oxygen, hydrogen, helium and nitrogen molecule (N2) 
%% densities, and temperature. The arrays produced contain each of the parameters
%% at 350km altitude, as well as the exospheric temperature and sigma shaping 
%% profile for the temperature profiles. The geographic arrays are assumed to be 
%% of the format where each row corresponds to the single UT in timmat at the 
%% same row index.
%% 
%% 9/3/2002 - add option to either call the MSIS-90 version, or the newer NRLMSIS2000.
%%

function [oxnt350,oxnt350c,hynt350,hynt350c,hent350,hent350c,n2nt350,n2nt350c,tmpn350,tmpn350c,tmpnex,tmpnexc,tmpnsp,tmpnspc]=modmsis3(dno,year,timmat,nopf,sms)

% dno = day number of interest (scalar)
% year = year of interest (scalar)
% timmat = UT time vector
% nopf = number of field lines to undertake full calculations
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% 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 = conjugate neutral hydrogen density at 350 km (m^-3)
% hent350 = neutral helium density at 350 km (m^-3)
% hent350c = conjugate neutral helium density at 350 km (m^-3)
% n2nt350 = neutral N2 density at 350 km (m^-3)
% n2nt350c = conjugate neutral N2 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 = conjugate neutral exospheric temperature (K)
% tmpnsp = neutral temperature shaping parameter
% tmpnspc = conjugate neutral temperature shaping parameter

re=6371;                                %Earth radi (km)
tp=length(timmat);                      %number of time points

%Calculate required solar parameters
[kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(1:3:22,dno*ones(1,8),year*ones(1,8));

date=(year-1900)*1e3+dno;               %required date format

if sms==0
 cd msis90
else
 cd msis2000
end

%Save parameters to disk
fid=fopen('data.txt','w');
fprintf(fid,'%5d %3d %4d %4d %5d\r\n',[date;tp;f107a(1);f107p(1);nopf]);
fprintf(fid,'%3d %3d %3d %3d %3d %3d %3d %3d\r\n',squeeze(apmat));
fclose(fid);

%call Fortran subrountine in MSIS directory to calculate required parameters
if sms==0
 !mlmsis93                    
else
 !mlmsis03                    
end

%load in results from MSIS
fid=fopen('res.txt','r');
resm=fscanf(fid,'%g %g %g %g %g %g %g',[7 inf]);
fclose(fid);

delete res.txt

cd ..

%record required results
pp1=1:(tp*nopf);
pp2=(tp*nopf+1):(2*tp*nopf);

tmpn350=resm(1,pp1)';
tmpnex=resm(2,pp1)';
tmpnsp=resm(3,pp1)';
oxnt350=1e6*resm(4,pp1)';
hynt350=1e6*resm(5,pp1)';
hent350=1e6*resm(6,pp1)';
n2nt350=1e6*resm(7,pp1)';

tmpn350c=resm(1,pp2)';
tmpnexc=resm(2,pp2)';
tmpnspc=resm(3,pp2)';
oxnt350c=1e6*resm(4,pp2)';
hynt350c=1e6*resm(5,pp2)';
hent350c=1e6*resm(6,pp2)';
n2nt350c=1e6*resm(7,pp2)';

tmpn350=reshape(tmpn350,nopf,tp);
tmpnex=reshape(tmpnex,nopf,tp);
tmpnsp=reshape(tmpnsp,nopf,tp);
oxnt350=reshape(oxnt350,nopf,tp);
hynt350=reshape(hynt350,nopf,tp);
hent350=reshape(hent350,nopf,tp);
n2nt350=reshape(n2nt350,nopf,tp);

tmpn350c=reshape(tmpn350c,nopf,tp);
tmpnexc=reshape(tmpnexc,nopf,tp);
tmpnspc=reshape(tmpnspc,nopf,tp);
oxnt350c=reshape(oxnt350c,nopf,tp);
hynt350c=reshape(hynt350c,nopf,tp);
hent350c=reshape(hent350c,nopf,tp);
n2nt350c=reshape(n2nt350c,nopf,tp);
