%% modiri5.m calls the IRI2000 to calculate electron density at the requested 
%% points of interest, as well as the HmF2. These points will be below the
%% boundary specified by the input. For points above the boundary,
%% 0 is returned.
%%
%% 31/7/2001: Updated to use IRI2001.
%%

function [eden,nmf2,hmf2,beta,eta,delta,zeta]=modiri5(lonmat,latmat,timmat,day,month,year,dismat,ionbd,inp,gms)

% lonmat = geographic longitude values matrix (rads)
% latmat = geographic latitude values matrix (rads)
% timmat = universal time matrix (fractional hours)
% day = day matrix
% month = month matrix
% year = year matrix
% dismat = displacement matrix, the distance from the centre of the Earth (km) 
% ionbd = upper altitude boundary on required ionosphere calculations (km) (scalar)   
% inp = IRI2000 FoF2 data set to be used (0 = CCIR, 1 = URSI) (scalar)
% gms = include IRI2001 geomagnetic storm model (0 = no, 1 = yes) (scalar)
%
% eden = electron density (e-/m^3)
% nmf2 = F2 electron density (e-/m^3)
% hmf2 = F2 height (km)
% beta, eta, delta, zeta = Bent model parameters

re=6371;                          %Earth radi (km)

%% Results matrix
eden=zeros(size(lonmat));
nmf2=eden;
hmf2=eden;
beta=eden;
eta=eden;
delta=eden;
zeta=eden;

%% First determine which points we need to use IRI2000 to calculate
ionc=find((dismat-re)<=ionbd);

if ~isempty(ionc)

 yeari=year(ionc);                %required year values
 dayi=day(ionc);                  %required day values
 monthi=month(ionc);              %required month values
 loni=lonmat(ionc);               %required longitude values
 lati=latmat(ionc);               %required latitude values
 timi=timmat(ionc);               %required time values
 alti=dismat(ionc)-re;            %required altitudes

 loni=longch2(loni);              %convert negative longitudes into positive values (add 2pi, required format of IRI2001)
 dnoi=daynum(dayi,monthi,yeari);  %determine day number 

 yeari=yeari(:);
 dnoi=dnoi(:);
 lati=lati(:);
 loni=loni(:);
 timi=timi(:);
 alti=alti(:);

 %Because of the delays in loading Ap values by IRI2000 when strom model used, sort by timmat before 
 %inputting so that points with same UT will be calculated sequentially. Note that this will only speed
 %up computational time if the day and year are the same for a given UT. Otherwise this sorting doesn't
 %actually help calculation time.
 [timi2,ti]=sort(timi);            %sort UT
 yeari2=yeari(ti);
 dnoi2=dnoi(ti);
 lati2=lati(ti);
 loni2=loni(ti);
 timi2=timi(ti);
 alti2=alti(ti);
 
 cd iri2000
 
 %For points for which the IRI2000 is required, save to file the required parameters
 fid=fopen('data.txt','w');
 fprintf(fid,'%d %d\r\n',[inp gms]);      
 fprintf(fid,'%4d %4d %4.3f %4.3f %6.4e %8.6e\r\n',[yeari2';dnoi2';lati2';loni2';timi2';alti2']);
 fclose(fid);
 
 %call Fortran subrountine in IRI2000 directory to calculate electron density
 !pwirid1                     

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

 delete res.txt 

 cd ..
 
end

edeni=zeros(1,size(resm,2));                  
nmf2i=edeni;                   
hmf2i=edeni;       
betai=edeni;
etai=edeni;
deltai=edeni;
zetai=edeni;

edeni(ti)=resm(1,:);                   
nmf2i(ti)=resm(2,:);                   
hmf2i(ti)=resm(3,:);                   
betai(ti)=resm(4,:);
etai(ti)=resm(5,:);
deltai(ti)=resm(6,:);
zetai(ti)=resm(7,:);

eden(ionc)=edeni;                  %required densities
nmf2(ionc)=nmf2i;                  %required NmF2 values
hmf2(ionc)=hmf2i;                  %required HmF2 values
beta(ionc)=betai;
eta(ionc)=etai;
delta(ionc)=deltai;
zeta(ionc)=zetai;

