%% iriap10.m is used to produce a look up table arrays of the electron 
%% density at 350km, HmF2, NmF2, HmF1 and the parameters for the Bent
%% model of the topside ionosphere for the given day of interest. The
%% arrays produced contain each of the parameters at 350km altitude.
%% 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.
%% 
%% 31/7/2001: Modified to call IRI2000.
%%

function [eden350,eden350c,hmf2,hmf2c,nmf2,nmf2c,beta,betac,eta,etac,delta,deltac,zeta,zetac]=iriap10(dno,year,timmat,nopf,inp,gms)

% dno = day number of interest (scalar)
% year = year of interest (scalar)
% timmat = UT time vector (fractional hours)
% nopf = number of field lines to undertake full calculations
% inp = FoF2 data set to be used (0 = CCIR, 1 = URSI) (scalar)
% gms = include IRI2000 geomagnetic storm model (0 = no, 1 = yes) (scalar)
%
% eden350 = electron density at 350km (m^-3)
% hmf2 = HmF2 (km)
% nmf2 = NmF2 (m^-3)
% beta, eta, delta, zeta = Bent model parameters
% #c = conjugate terms

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

cd iri2000

%Save parameters to disk
fid=fopen('data.txt','w');
fprintf(fid,'%4d %3d %3d %5d %d %d\r\n',[year;dno;tp;nopf;inp;gms]);
fclose(fid);

%call Fortran subrountine in IRI2000 directory to calculate electron density
!pwirid3                   
 
%load in results
fid=fopen('res.txt','r');
resm=fscanf(fid,'%e %g %e %e %e %e %e',[7 inf]);
fclose(fid);

delete res.txt

cd ..

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

eden350=resm(1,pp1)';
hmf2=resm(2,pp1)';
nmf2=resm(3,pp1)';
beta=resm(4,pp1)';
eta=resm(5,pp1)';
delta=resm(6,pp1)';
zeta=resm(7,pp1)';

eden350c=resm(1,pp2)';
hmf2c=resm(2,pp2)';
nmf2c=resm(3,pp2)';
betac=resm(4,pp2)';
etac=resm(5,pp2)';
deltac=resm(6,pp2)';
zetac=resm(7,pp2)';

eden350=reshape(eden350,nopf,tp);
hmf2=reshape(hmf2,nopf,tp);
nmf2=reshape(nmf2,nopf,tp);
beta=reshape(beta,nopf,tp);
eta=reshape(eta,nopf,tp);
delta=reshape(delta,nopf,tp);
zeta=reshape(zeta,nopf,tp);

eden350c=reshape(eden350c,nopf,tp);
hmf2c=reshape(hmf2c,nopf,tp);
nmf2c=reshape(nmf2c,nopf,tp);
betac=reshape(betac,nopf,tp);
etac=reshape(etac,nopf,tp);
deltac=reshape(deltac,nopf,tp);
zetac=reshape(zetac,nopf,tp);
