% This routine searches many time series of monthly climate data
% for changes to the intensity and period of annual cycle.
% Output can be stored and sorted to find the time series (land 
% cell) that shows the most variance in its annual cycles.
% Adjusting the frequency range (f) expands or contracts the range 
% to be analyzed. However, there is a performance trade-off, as 
% broader ranges require more comutation.
%
close all;clear
l=2^(16);% Transform length. FFT will be fastest when this is a power of 2.
         % This could be decreased to speed up compuatation at the expense
         % of accuracy of the peak's location.
nfields=3;
nelements=120;% The number of elements or elements per decade
i=1;
% search serially through grid cell output files
%for gc=0:2795
for gc=1:1
    fname=sprintf('pcm_B05_12_hpt_%04i_monthly.dat', gc);
%    printf('\n');
%    printf('Loading file: %s\n', fname);
    a(:,:,i)=load(fname);
    for j=1:nfields
        s(:,1)=a(:,j:j,i);
        f=[0:1:round(0.1*l)]/l;
        for k=0:8
            s(1:nelements,k+1)=s(k*nelements+1:k*nelements+nelements,1);% Window
            sd_s=std(s(1:nelements,k+1));% Normalize to variance of 1
            s_n=s(1:nelements,k+1)/sd_s;
            s_ft=fft(s_n,l);% Do FFT
            sp=s_ft.*conj(s_ft)/l;% Power of the DFT
            [a_k cnt]=max(sp(4501:6000));% Declare variables for real (a_k)
                                         % and index (cnt)
%            fprintf('Peak amp: %0.5f, Period (months): %0.2f\n', a_k, 1/f(4500+cnt))
            pa(k+1,j)=a_k;
            pe(k+1,j)=1/f(4500+cnt);
        end
    end
%    fprintf('Peak amplitude diff 1st/last decades: %0.5f, Variance: %0.5f\n', pa(1,1)-pa(9,1), var(pa,1)) %change in peak amplitude
%    Output all peak amplitudes & their periods to array gpa_pe
    gc_cnt=gc*9;
    gpa_pe(gc_cnt+1:gc_cnt+9,1)=gc+1;
    gpa_pe(gc_cnt+1:gc_cnt+9,2)=pa(:,1);
    gpa_pe(gc_cnt+1:gc_cnt+9,3)=pa(:,2);
    gpa_pe(gc_cnt+1:gc_cnt+9,4)=pa(:,3);
    gpa_pe(gc_cnt+1:gc_cnt+9,5)=pe(:,1);
    gpa_pe(gc_cnt+1:gc_cnt+9,6)=pe(:,2);
    gpa_pe(gc_cnt+1:gc_cnt+9,7)=pe(:,3);
    % Global difference (first - last) in amplitude for all 3 vars
    gda(gc+1,1)=(pa(1,1)-pa(9,1))/pa(1,1)/.01;
    gda(gc+1,2)=(pa(1,2)-pa(9,2))/pa(1,2)/.01;
    gda(gc+1,3)=(pa(1,3)-pa(9,3))/pa(1,3)/.01;
    % Normalized (to variance of 1) global variance in amplitude for all 3 vars
    gva(gc+1,1)=var(pa(:,1))/std(pa(:,1));
    gva(gc+1,2)=var(pa(:,2))/std(pa(:,2));
    gva(gc+1,3)=var(pa(:,3))/std(pa(:,3));
    % fprintf('Period diff 1st/last decades %0.5f, Variance %0.5f\n', pe(1,1)-pe(9,1), var(pe,1)) %change in period
    % Global difference (first - last) in period for all 3 vars.
    % Signs are switched so that an increase in period length between
    % the first and last decade will show up as a positive number.
    gdp(gc+1,1)=-(pe(1,1)-pe(9,1))/pe(1,1)/.01;
    gdp(gc+1,2)=-(pe(1,2)-pe(9,2))/pe(1,2)/.01;
    gdp(gc+1,3)=-(pe(1,1)-pe(9,3))/pe(1,3)/.01;
    % Normalized (to variance of 1) global variance in period for all 3 vars
    gvp(gc+1,1)=var(pe(:,1))/std(pe(:,1));
    gvp(gc+1,2)=var(pe(:,2))/std(pe(:,2));
    gvp(gc+1,3)=var(pe(:,3))/std(pe(:,3));
    ++i;
end
fprintf('Global Difference (first decade - last decade) in Amplitude as a percentage of the first decades''s amplitude. Negative percentages indicate reduced intensity.\n')
gda
fprintf('Global Variance in Amplitude as a percentage of the maximum amplitude.\n')
gva
fprintf('Global Difference (first decade - last decade) in Period as a percentage of the first decades''s period. Negative percentages indicate shorter period.\n')
gdp
fprintf('Global Variance in Period as a percentage of the maximum period.\n')
gvp
%
fprintf('Max/Min Global Difference in Amplitude. Subtract 1 from grid cell number to identify grid cell data file.\n')
fprintf('Variable 1\n')
[gdamax grid_cell]=max(gda(:,1))
[gdamin grid_cell]=min(gda(:,1))
fprintf('\n')
fprintf('Variable 2\n')
[gdamax grid_cell]=max(gda(:,2))
[gdamin grid_cell]=min(gda(:,2))
fprintf('\n')
fprintf('Variable 3\n')
[gdamax grid_cell]=max(gda(:,3))
[gdamin grid_cell]=min(gda(:,3))

save('pcm_locate.out', 'gpa_pe', 'gda', 'gva', 'gdp', 'gvp')
