% The purpose of this code is to animate in time the work of Robert M
% Nadeau on the repeating earthquake sequences of the creeping region of
% the central San Andreas Fault.  He has created a figure which has been
% updated and used in various forms over the last several years to
% illustrate the "pulsing" of the creep rate in the area.  The figure in
% question can be found in the paper "Periodic Pulsing of Characteristic
% Microearthquakes on the San Andreas Fault" by Nadeau and McEvilly, 2004.
% (http://www.sciencemag.org/content/303/5655/220.full)  The figure we are
% attempting to animate is the left pane of figure two.  
%
% It is our hope to make the data in the figure more immediately intuitive
% so that the observer may get a more direct impression of the slip
% transients in the region and more readily comprehend the data.
%
% Created by Ryan Turner for a UC Berkeley class on MATLAB programming,
% EPS109.  
%

% Clear variables and figures
clear;
clf;

% Import tabled data used to create the figure as a matrix
load DB.stp2.HDD;
A=DB_stp2(:,2);     % Pulls out the Y-coordinate (along strike)
B=DB_stp2(:,4);     % Pulls out the decimal date of event
C=DB_stp2(:,7);     % Pulls out the incremental slip associated with event
D=horzcat(B,A,C);    % Concatenates the data into a (Y,DATE,ISLIP) matrix

% Sort the matrix along the date column (1)
S=sortrows(D);
Ne = size(D,1);

% Set Variables
dt = 0.01;     %One decimal day is ~0.0027397
T = round(max(B)-min(B))+1;     %About 27 years of data
Nt = round(T/dt);

%Determine number of sites and index
Au = sort(unique(A));
Ny = size(Au,1);
P = zeros(Ny,1);

%Set initial time
ti=1984;
tb = 1984:dt:2012;
tbh = zeros(size(tb,2),1);
f = 1;  %start the frame counter

%Creates a matrix of time bins associated with each of event entry in S
for h = 1:Ne;   %for each event entry (~3500)
    t=ti;       %restart clock
    for i = 1:Nt    %for each time bin
        t2 = t+i*dt;    %set end time of bin
        if S(h,1) >= t && S(h,1) < t2   %if entry has date >= t, <t2
            break   %break out of for loop to ***
        end
    end    
    tbh(h) = tb(i); %***assign the matching bin to that entry
    t = t2;         %move to next time bin and resume for loop over i
end    

for i = 1:size(tb,2);   %for each time bin...
    for h = 1:Ne;   %and for each entry (~3500)
        if tbh(h) >= tb(i) && tbh(h) < tb(i+1); %check if entry is in time bin
            for j = 1:Ny-1;   %if entry is in time bin, see which site it goes to
                if S(h,2) >= Au(j) && S(h,2) < Au(j+1); %when site is found
                    P(j) = P(j) + S(h,3);   %write P value for that time and site
                    break
                end
            end
        end
    end
    t = tb(i);
    %plot(P,Au,'k-')
    area(P,Au,'FaceColor',[.5 .5 .5],'EdgeColor',[0 0 0],'LineWidth',2);
    colormap gray;
    title([sprintf('Date = %4.3f',t)]);
    axis([0 160 min(Au) max(Au)]);
    xlabel('Cumulative Slip [cm]')
    ylabel('Distance Along Fault Strike From Middle Mountain')
    drawnow
    %pause(0.1)
    saveas(gcf, ['C:\Temp\movie\slip_animation_',sprintf('%05d',f),'.png'])
    f = f+1;    %advance frame number
end

