%% iriap8.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.
%% iriap8.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. 
%%
%% 31/7/2001: Modified to call IRI2000.
%% 

function [eden350,eden350c,hmf2,hmf2c,nmf2,nmf2c,beta,betac,eta,etac,delta,deltac,zeta,zetac]=iriap8(lon350,lat350,lon350c,lat350c,dno,year,dtv,inp,gms)

% 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)
% inp = IRI2000 FoF2 data set to be used (0 = CCIR, 1 = URSI) (scalar)
% gms = include IRI2001 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)

%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

lon350=longch2(lon350);                 %convert negative longitudes into positive values (add 2pi, required format of IRI2001)
lon350c=longch2(lon350c);               %convert negative longitudes into positive values (add 2pi, required format of IRI2001)

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

cd iri2000

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

%Save geographic co-ordinates to disk
fid=fopen('data2.txt','w');
fprintf(fid,'%4.3f %4.3f\r\n',[[lon350(:)' lon350c(:)'];[lat350(:)' lat350c(:)']]);
fclose(fid);

%call Fortran subrountine in IRI2000 directory to calculate electron density
!pwirid2                   

%load in results from IRI2000
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
resm2=reshape(resm,(sch*7),(tp*2));

resm3=resm2(:,(1:2:(tp*2)));
resm3c=resm2(:,(2:2:(tp*2)));

pp1=(1:7:(sch*7))';

eden350=resm3(pp1,:);
hmf2=resm3((pp1+1),:);
nmf2=resm3((pp1+2),:);
beta=resm3((pp1+3),:);
eta=resm3((pp1+4),:);
delta=resm3((pp1+5),:);
zeta=resm3((pp1+6),:);

eden350c=resm3c(pp1,:);
hmf2c=resm3c((pp1+1),:);
nmf2c=resm3c((pp1+2),:);
betac=resm3c((pp1+3),:);
etac=resm3c((pp1+4),:);
deltac=resm3c((pp1+5),:);
zetac=resm3c((pp1+6),:);
