%% f107day.m determines the daily F10.7cm flux value for the current day, the
%% previous day and the average from 81 days centered on the required day 
%% (required by MSIS-86). This is done from the data downloaded from the NATIONAL 
%% GEOPHYSICAL DATA CENTER.
%% The day value is a 'running mean' based on the next day and the the current day 
%% (if UT>12) or the current day and previous day (UT<12). The previous day F107
%% is calculated in a similar manner.
%% 

function [f107,f107p,f107av]=f107day(timmat,dno,year)

% timmat = UT matrix
% dno = day number matrix
% year = year matrix
% f107 = daily solar flux
% f107p = daily solar flux - previous day
% f107av = 81 day average solar flux

sz=size(dno);                           %size of input
f107=zeros(sz);                         %f107 results matrix
f107p=f107;                             %previous day results matrix
f107av=f107;                            %f107 average results matrix
nop=sz(1)*sz(2);                        %number of points of interest

%convert input to column vector
timmat=timmat(:);
dno=dno(:);
year=year(:);

%Previous year
yearo=-1;

for i=1:nop
 timi=timmat(i);                                  %i'th UT
 dnoi=dno(i);                                     %i'th day number
 yeari=year(i);                                   %i'th year
 if yeari~=yearo
  eval(['load data\f107day\f',int2str(yeari)])    %load in data for i'th year
  eval(['f7d=f',int2str(yeari),';'])              %relabel i'th data f7d
  eval(['clear f',int2str(yeari)])                %clear i'th data 
  yearo=yeari;                                    %record new previous year value
 end 
 
 %Determine required i'th F10.7
 if timi==12
  f107(i)=f7d(dnoi,4);                             
 elseif timi<12
  m=(f7d(dnoi,4)-f7d(dnoi,3))/24;
  f107(i)=f7d(dnoi,3)+m*(timi+12);
 else
  m=(f7d(dnoi,5)-f7d(dnoi,4))/24;
  f107(i)=f7d(dnoi,4)+m*(timi-12);
 end 

 %Determine required i'th previous F10.7
 if timi==12
  f107p(i)=f7d(dnoi,3);                             
 elseif timi<12
  m=(f7d(dnoi,3)-f7d(dnoi,2))/24;
  f107p(i)=f7d(dnoi,2)+m*(timi+12);
 else
  m=(f7d(dnoi,4)-f7d(dnoi,3))/24;
  f107p(i)=f7d(dnoi,3)+m*(timi-12);
 end 
 
 f107av(i)=f7d(dnoi,6);                           %required i'th 81 day average F10.7
end 

