% Simple demo code to illustrate process of loading a matlab binary file 
% and processing a set seismograms.
%
% Set your MATLAB working directory to a directory containing SAC files.
% This code will read in the SAC files, and tape, bandpass, and decimate
% the files.
%
% Created for the USArray Data Processing and Analysis Short Course
% August 15-19, 2011, Northwester University
% author: % Bob Woodward, IRIS
%
clear
%
tic
load mydata
toc
%
% preallocate memory for speed
[m n] = size(sacdat);
envdat = zeros(m,1);
tdiff  = zeros(n,1);
%
%
% Processing loop
%
% We pick time series #1 as our reference and cross correlate it against
% every other time series in the matrix.
%
env_ref = hilbert2(sacdat(1:npts(1),1),1/dt(1));

for i=1:length(lat)
    env = hilbert2(sacdat(1:npts(i),i),1/dt(i));
    np = min([npts(1) npts(i)]);
    cc = xcorr(env_ref(1:np), env(1:np) ); % compute the cross correlation
    [x n] = max(abs(cc));   % extract peak value and array index of peak value
    tdiff(i) = tbeg(1) - tbeg(i) - ( n - np ) * dt(1); % compute time shift
end;
%
% plot some results
%
hold on
% load and plot coastline and political boundaries
shore     = importdata('shoreline_US.dat');
political = importdata('political_US.dat');
plot(shore(:,1),shore(:,2),'k');
plot(political(:,1),political(:,2),'k');
%
% plot time offset to peak correlation
set(gca,'CLim',[-800 0]); % scale for scatter plot
scatter(lon,lat,100,tdiff,'filled','MarkerEdgeColor','k');
% put a red ring around our reference station
scatter(lon(1),lat(1),300,'MarkerEdgeColor','r');
colorbar % add a colorbar to plot
axis tight

