%% apmsis4.m is a modified version of apmsis3.m. 
%% As with apmsis3.m it calculates the current 3-hourly Ap, the mean of the 
%% current plus seven previous 3-hourly Ap and the minimum of the current and 
%% next seven 3-hourly Ap. Further, it now calculates the Ap values required
%% by the MSIS model which 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
%% Note that the daily Ap will be taken to be equal to the mean of the current and
%% seven previous 3-hourly Ap values (ie the 2nd Ap calculated). The current 3-hourly
%% Ap is already outputted (the 1st Ap calculated). 
%% Each row is for a 3 hour period, with the columns being day number,
%% UT, and the Ap values. The files are saved every year, with the 
%% naming format being apm# where # = year. Thus the outputted Ap values are:
%% 1 - current 3-houlry Ap.
%% 2 - mean of the current plus seven previous 3-hourly Ap.
%% 3 - minimum of the current and next seven 3-hourly Ap.
%% 4 - (3) above
%% 5 - (4) above
%% 6 - (5) above
%% 7 - (6) above
%% 8 - (7) above
%% 

clear

years=input('Starting File Year? ');
yearf=input('Finishing File Year? ');

nop=yearf-years+1;                           %number of files to be converted
ap27=zeros(1,27);                            %next 7 3-hourly Ap, current Ap and last 19 3-hour Ap values

for i=1:nop
 if i==1
  yeari=years;                          %year index
 else
  yeari=yeari+1;                        %increment year index
 end 
 
 eval(['load a',int2str(yeari)])        %load in i'th year data file
 ap3h=ap3h';                            %take transpose of 3-hour Ap matrix data
  
 dn=0;                                  %current day number
 h3i=0;                                 %current 3-hour index
 jmax=8*length(apv);                    %number of 3-hour periods in current year
 apmat=zeros(jmax,10);                  %results matrix for i'th year

 for j=1:jmax                           %loop through each day of i'th year
  if rem(j-1,8)==0                      %check if first 3-hour period of new day
   dn=dn+1;                             %update day number
   timut=0;                             %reset UT
  else 
   timut=timut+3;                       %update new UT
  end
   
  ap27(2:27)=ap27(1:26);                %move 3-hour Ap values back one time interval
  h3i=h3i+1;                            %increment h3i by one
  ap27(1)=ap3h(h3i);                    %j'th 3-hour Ap increment
  
  %Record the j'th Ap values
  apmat(j,:)=[dn timut ap27(8) mean(ap27(8:15)) min(ap27(1:8)) ap27(9) ap27(10) ap27(11) mean(ap27(12:19)) mean(ap27(20:27))];   
  
 end
 
 %Save the results for the i'th year
 eval(['apm',int2str(yeari),'=apmat;'])
 eval(['save apm',int2str(yeari),' apm',int2str(yeari)])  
 eval(['clear apm',int2str(yeari)]) 
end 

%% Now need to loop through and correct data which is offset by one day
for i=1:nop
 if i~=1                                %update year if necessary
  yeari=yeari-1;
 end
 
 %Load the results for the i'th year
 eval(['load apm',int2str(yeari),' apm',int2str(yeari)])  
 eval(['nopd=size(apm',int2str(yeari),',1)/8;'])
 
 %Remove first day of data, which belongs as last day from the previous year
 eval(['apmco=apm',int2str(yeari),'(1:8,:);'])
 eval(['apm',int2str(yeari),'(1:8,:)=[];'])
 
 %Update the day number
 eval(['apm',int2str(yeari),'(:,1)=apm',int2str(yeari),'(:,1)-1;'])
 
 %Update the last day with its correct value
 if i==1
  apmco2=[nopd*ones(8,1) zeros(8,9)];
  eval(['apm',int2str(yeari),'=[apm',int2str(yeari),';apmco2];'])
 else
  apmco2(:,1)=nopd*ones(8,1);
  eval(['apm',int2str(yeari),'=[apm',int2str(yeari),';apmco2];'])
 end 
 
 %Update ampco
 apmco2=apmco;
 
 %Save modified files
 eval(['save apm',int2str(yeari),' apm',int2str(yeari)])  
 eval(['clear apm',int2str(yeari)])   
end 