%% msisap9.m is calls the MSIS-90 model to calculate neutral parameters oxygen
%% density, hydrogen density 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.
%% 

function [oxnt350,oxnt350c,hynt350,hynt350c,tmpn350,tmpn350c,tmpnex,tmpnexc,tmpnsp,tmpnspc]=msisap9(dno,year,timmat,nopf)

% dno = day number of interest (scalar)
% year = year of interest (scalar)
% timmat = UT time vector
% nopf = number of field lines to undertake full calculations
%
% 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

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

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

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

cd msis90

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

%call Fortran subrountine in MSIS90 directory to calculate required parameters
!mlmsis93                    

%load in results from MSIS90
fid=fopen('res.txt','r');
resm=fscanf(fid,'%g %g %g %g %g',[5 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)';

tmpn350c=resm(1,pp2)';
tmpnexc=resm(2,pp2)';
tmpnspc=resm(3,pp2)';
oxnt350c=1e6*resm(4,pp2)';
hynt350c=1e6*resm(5,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);

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