%% msisap7.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. msisap7.m expects 
%% the inputs for the longitude and latitudes to be column vectors, and converts
%% to such if they are not. The output arrays have size length(lon350) by dtv. 
%% 

function [oxnt350,oxnt350c,hynt350,hynt350c,tmpn350,tmpn350c,tmpnex,tmpnexc,tmpnsp,tmpnspc]=msisap7(lon350,lat350,lon350c,lat350c,dno,year,dtv)

% lon350 = geographic longitude at 350km (rads)
% lat350 = geographic latitude at 350km (rads)
% lon350c = conjugate geographic longitude at 350km (rads)
% lat350c = conjugate geographic latitude at 350km (rads)
% dno = day number of interest (scalar)
% year = year of interest (scalar)
% dtv = UT time step to use which must divide exactly into 24 (fractional hours,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 = 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)

%convert input to column matrices
lon350=lon350(:);
lat350=lat350(:);
lon350c=lon350c(:);
lat350c=lat350c(:);

tp=round(24/dtv);                       %number of time points
sch=length(lon350);                     %length of lon350
alt350=350;                             %350km altitude   

%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

%convert from rads to degrees
lon350=rad2deg(lon350);                 
lat350=rad2deg(lat350);                 
lon350c=rad2deg(lon350c);                 
lat350c=rad2deg(lat350c);                 

cd msis90

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

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

%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*sch);
pp2=(tp*sch+1):(2*tp*sch);

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,tp,sch)';
tmpnex=reshape(tmpnex,tp,sch)';
tmpnsp=reshape(tmpnsp,tp,sch)';
oxnt350=reshape(oxnt350,tp,sch)';
hynt350=reshape(hynt350,tp,sch)';

tmpn350c=reshape(tmpn350c,tp,sch)';
tmpnexc=reshape(tmpnexc,tp,sch)';
tmpnspc=reshape(tmpnspc,tp,sch)';
oxnt350c=reshape(oxnt350c,tp,sch)';
hynt350c=reshape(hynt350c,tp,sch)';
