%% ionfreq.m calculates the total ionisation frequency of either He+ or O+, up to 
%% the maximum wavelength at which He+ is produced. It also returns the individual
%% ionisation frequencies and cross sections for each bin. This done using the 
%% results from the Torr et al. GRL 1979 results and the EUVAC model.
%%

function [qtot,qi,cs]=ionfreq(f107,f107a,ionr)

% f107 = daily F10.7 of day of interest (scalar)
% f107a = 81 day averaged F10.7 centered on the day of interest (scalar)
% ionr = ion required (0 = He+, 1 = O+)
%
% qtot = total ionisation frequency (s^-1)
% qi = ionisation frequency for each bin (s^-1)
% cs = cross section of each bin (cm^2)

%Total O+ cross section (cm^2)
ocs=1e-18*[1.06 3.53 5.96 7.55 8.43 9.26 8.78 9.70 9.72 10.03 10.84 10.70 11.21 11.25 11.64 11.91];

%Total He+ cross section  (cm^2)
hecs=1e-18*[0.21 0.53 1.02 1.71 2.16 2.67 2.38 3.05 3.05 3.65 4.35 4.25 5.51 6.53 7.09 0.72];

%Calculate the solar flux in each bin
solf=euvac(f107,f107a);

%Total ionisation frequency is given by the sum of the solar flux multiplied by the cross section.
if ionr==0
 qi=hecs.*solf;
 qtot=sum(qi);
 cs=hecs; 
else
 qi=ocs.*solf;
 qtot=sum(qi);
 cs=ocs;
end