% 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, = 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