%% modmsis1.m calculates the neutral oxygen, hydrogen, helium and nitrogen
%% molecule (N2) densities, the neutral temperature at the required altitude, 
%% the exospheric temperature and the temperature shaping parameter using the
%% neutral atmospheric MSIS model. 
%% modmsis1.m is based on the earlier modmsis9.m.
%%
%% 27/10/2001 - Modified to use 7 element Ap array input. Note that if the
%% inputs are size m x n, then apday needs to be m x n x 7, which is the
%% output size given by solarval.m
%% 
%% 12/11/2001 - modify so that if results required at one or more locations, but 
%% all at the same UT and date, then user can input timmat, day, month, year, f107a,
%% and f107p as scalars, and apmat as a single 1x1x7 vector.
%% 
%% 9/3/2002 - add option to either call the MSIS-90 version, or the newer NRLMSIS2000.
%%

function [oxnt,hynt,hent,n2nt,tmpn,tmpnex,tmpnsp]=modmsis1(lonmat,latmat,dismat,timmat,day,month,year,f107a,f107p,apmat,sms)

% lonmat = geographic longitude values matrix (rads)
% latmat = geographic latitude values matrix (rads)
% dismat = displacement matrix, the distance from the centre of the Earth (km)
% timmat = universal time (UT) matrix (fractional hours) (vector or scalar - see 12/11/2001 note)
% day = day matrix (vector or scalar - see 12/11/2001 note)
% month = month matrix (vector or scalar - see 12/11/2001 note)
% year = year matrix (vector or scalar - see 12/11/2001 note)
% f107a = average f107 3 month value (vector or scalar - see 12/11/2001 note)
% f107p = daily f107 value from previous day (vector or scalar - see 12/11/2001 note)
% apmat = 7 Ap values required by MSIS (see note 27/10/2001) (vector or scalar - see 12/11/2001 note)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% oxnt = neutral oxygen density (m-3)
% hynt = neutral hydrogen density (m-3)
% hent = neutral helium density (m-3)
% n2nt = neutral N2 density (m-3)
% tmpn = neutral temperature (K)
% tmpnex = neutral exospheric temperature (K)
% tmpnsp = neutral temperature shaping parameter

re=6371;                                %Earth radi (km) 
sch=size(dismat);                       %input size

dno=daynum(day,month,year);             %calculate day number
date=(year-1900)*1e3+dno;               %required date format
timl=timut2l(timmat,lonmat);            %calculate local time 
lonmat=rad2deg(lonmat);                 %convert from rads to degrees
latmat=rad2deg(latmat);                 %convert from rads to degrees
timmat=3600*timmat;                     %convert to seconds
dismat=dismat-re;                       %convert to altitude        

date=date(:)';
timmat=timmat(:)';
dismat=dismat(:)';
latmat=latmat(:)';
lonmat=lonmat(:)';
timl=timl(:)';
f107a=f107a(:)';
f107p=f107p(:)';
if prod(size(date))==1
 apmat2=squeeze(apmat);
else
 apmat2=reshape(apmat,prod(sch),7);
 apmat2=apmat2';
end

if sms==0
 cd msis90
else
 cd msis2000
end

%Save parameters to disk
if prod(size(date))==1
 fid=fopen('data.txt','w');                           
 fprintf(fid,'%5d %5d %4d %4d %3d %3d %3d %3d %3d %3d %3d\r\n',[date;timmat;f107a;f107p;apmat2]);
 fprintf(fid,'%8.6e %4d %4d %6.4e\r\n',[dismat;latmat;lonmat;timl]);
 fclose(fid);
 %call Fortran subrountine in MSIS directory to calculate required parameters
 if sms==0
  !mlmsis91b                    
 else
  !mlmsis01b
 end
else
 fid=fopen('data.txt','w');                           
 fprintf(fid,'%5d %5d %8.6e %4d %4d %6.4e %4d %4d %3d %3d %3d %3d %3d %3d %3d\r\n',[date;timmat;dismat;latmat;lonmat;timl;f107a;f107p;apmat2]);
 fclose(fid);
 %call Fortran subrountine in MSIS directory to calculate required parameters
 if sms==0
  !mlmsis91                    
 else
  !mlmsis01
 end
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 ..

oxnt=1e6*reshape(resm(1,:),sch);
hynt=1e6*reshape(resm(2,:),sch);
hent=1e6*reshape(resm(3,:),sch);
n2nt=1e6*reshape(resm(4,:),sch);
tmpn=reshape(resm(5,:),sch);
tmpnex=reshape(resm(6,:),sch);
tmpnsp=reshape(resm(7,:),sch);




