%% solarval.m calculates the following solar parameters for the
%% the points of interest.
%% 1) kp = 3-hourly Kp solar activity
%% 2) kpmin = minimum 3-hourly Kp in the next 24 hours
%% 3) apday = daily Ap (average of latest eight 3-hourly Ap values)
%% 4) f107 = current F10.7 flux 
%% 5) f107p = F10.7 flux 24 hours previously 
%% 6) f107a = average 81 day F10.7 flux centered on current day
%% 7) apmat = 3-D matrix which contains the Ap values required by MSIS model. If input is
%%    size m x n, then apmat is m x n x 7, where the 7 values are:
%%   (1) DAILY AP
%%   (2) 3 HR AP INDEX FOR CURRENT TIME
%%   (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
%%   (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
%%   (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
%%   (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR TO CURRENT TIME
%%   (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 59 HRS PRIOR TO CURRENT TIME
%%      
%% Ap data is downloaded from ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/
%% and F10.7 is from NOAA /STP/SOLAR_DATA/SOLAR_RADIO/FLUX DAILYPLT.OBS
%%

function [kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timmat,dno,year)

% timmat = universal time (UT) matrix (fractional hours)
% day = day number matrix
% year = year matrix
%
% kp = 3-hourly Kp daily solar activity value
% kpmin = minimum 3-hourly Kp in the next 24 hours
% apday = 'running' daily Ap (mean of latest eight 3-hour periods)
% f107 = 'running ' daily f107 value (linear interpolation)
% f107p = 'running ' daily f107 value (linear interpolation of previous days value)
% f107a = average f107 3 month value
% apmat = Ap parameters required by MSIS

%Warning about changed solarval output
%disp(' ')
%disp('Warning: using solarval - have you included changed output format?')
%disp('Have a nice day :-)')
%disp(' ')

%Calculate required F10.7 values
[f107,f107p,f107a]=f107day(timmat,dno,year);

kp=zeros(size(timmat));                 %Kp results matrix
kpmin=kp;                               %Minimum Kp results matrix
apday=kp;                               %Ap results matrix
nop=prod(size(kp));                     %number of inputs
apmat=zeros([size(timmat) 7]);          %MSIS Ap results matrix

%Calculate the required Ap values
yearo=-1;                               %previous year
for i=1:nop
 timi=timmat(i);                        %i'th timmat                             
 dnoi=dno(i);                           %i'th day number
 yeari=year(i);                         %i'th year
 if yeari~=yearo                                       %check if i'th year from previous year
  eval(['load data\ap3hour\apm',int2str(yeari)])       %load in data for i'th year and month
  eval(['apd=apm',int2str(yeari),';'])                 %relabel i'th data apd
  eval(['clear apm',int2str(yeari)])                   %clear i'th data 
  yearo=yeari;                                         %update yearo
 end
 idh=floor(timi/3)+1;                                  %index of i'th time
 idh(idh>8)=8;                                         %check for 24hour value
 id=8*(dnoi-1)+idh;                                    %index of i'th Ap value
 apday(i)=apd(id,4);                                   %required i'th Ap value
 kp(i)=apd(id,3);                                      %required i'th Kp as Ap (convert below)
 kpmin(i)=apd(id,5);                                   %required i'th Minimum Kp as Ap (convert below)
 [ii,ij]=ind2sub(size(timmat),i);                      %row and column index of i'th point
 apmat(ii,ij,:)=[apd(id,4) apd(id,3) apd(id,6:10)];    %required i'th Ap values for MSIS 
end 

%Calculate Kp
kp=ap2kp(kp);                                          %convert latest 3-hourly Ap to 3-hourly Kp
kpmin=ap2kp(kpmin);                                    %convert minimum 3-hourly Ap to 3-hourly Kp