% Test program for animating a wavefield
%
% This version is a variant of the usual red and blue blinking dots. This
% version contours the wavefield - whihc is more useful for examining
% subtle features of the wavefronts.
%
% Bob Woodward
% Incorporated Research Institutions for Seismology / IRIS
% woodward@iris.edu
%
% This code illustrates the basics of one way to do a wavefield animation
% within MATLAB.
%
% Inspired by, and modeled after, the work of Chuck Ammon, Pennsylvania
% State University.
%
% Assumes data are in SAC files that are in MATLAB path.
% Works for BH or LH data - or other sample rates.
% Should work even if sample rates are mixed.
%
clear
%
% define colormap - load into array cmap -  just simple 64x3 array that
% expresses the pretty blue and red colors in r g b values.
%
% don't need dotcolors in this version - we aren't plotting dots
%load dotcolors;
%
% Read shoreline and political boundary data for map.
% You can get coastlines and political boundaries from NOAA at
% http://www.ngdc.noaa.gov/mgg_coastline/
% Use the pull down menu labeled "coastline database" to select either a
% coastline or political boundary database. Note the option for matlab
% format output.
%
shore     = importdata('shoreline_western_US.dat');
political = importdata('political_western_US.dat');
%
% To pre-process the seismograms with SAC, use commands like the following:
% (assuming the SAC filenames start with 2008*)
% After filtering I write all the files, but changing the prefix.
%
% SAC> r 2008*.SAC
% SAC> rmean
% SAC> taper
% SAC> lp co 0.02 np 4 be p 2
% SAC> w change 2008 test
%
% read directory
dirlist = dir('*.SAC');
fprintf('Number of SAC files: %d\n',length(dirlist));
% read file headers only
for i=1:length(dirlist)
    s = rsac_h(dirlist(i).name);
    maxpts = max(0,s.npts);
end;
fprintf('scanned file headers\n');
% preallocate for speed
sacdata = zeros(maxpts,length(dirlist));
lat  = zeros(length(dirlist),1);
lon  = zeros(length(dirlist),1);
tbeg = zeros(length(dirlist),1);
npts = zeros(length(dirlist),1);
dt   = zeros(length(dirlist),1);
amp  = zeros(length(dirlist),1);
%
% Read the data, extrac necessary metadata, filter, then store the 
% seismograms in the columns of matrix sacdata.
% Get station locations and start time of seismograms from SAC header.
% tbeg equals UNIX epochal time, in seconds.
% Demean and then filter the data.
% Don't use MATLAB array operation to compute mean, as array has a lot of 
% zero padding because the seismograms are of different length.
% We assume that s.scale is in units of counts / (m/s)
%
% The scaling is done in case we want to use absolute scaling across the
% network.
%
% Should really include tapering prior to filtering!
%
% The epochal time functions (epoch, ep2dat) could use MATLAB built-in
% functions instead (with a little work).
%
for i=1:length(dirlist)
    s = rsac(dirlist(i).name);
    lat(i) = s.stla;
    lon(i) = s.stlo;
    % compute begin time of seismogram
    tbeg(i) = 86400*datenum(s.nzyear, 0, s.nzjday, s.nzhour, s.nzmin, s.nzsec+s.nzmsec/1000.) + s.b;
    npts(i) = s.npts;
    dt(i)   = s.delta;
    % demean and convert to nm/s
    s.d(1:npts(i)) = s.d(1:npts(i)) - mean(s.d(1:npts(i)));
    s.d(1:npts(i)) = s.d(1:npts(i)) / s.scale * 1.e9; % converts to units of nm/s
    % now filter
    nq = 0.5 / dt(i); % Nyquist
    % high pass?
    co    = 0.001 / nq; % normalized freq = corner freq [Hz] / Nyquist [Hz]
    [b a] = butter(4,co,'high');
    s.d(1:npts(i)) = filter(b,a,s.d(1:npts(i)));
    % low pass?
    co    = .01 / nq; % normalized freq = corner freq [Hz] / Nyquist [Hz]
    [b a] = butter(4,co,'low');
    s.d(1:npts(i)) = filter(b,a,s.d(1:npts(i)));
    % put data in big array
    sacdata(1:npts(i),i) = s.d(1:npts(i));
end;
fprintf('read data and filtered\n');
%
% Set movie start time, duration, and time step
%
tstart = min(tbeg) + 0.0;
tend   = tstart + 2000.0;  % set movie duration, in s
tstep  = 10.0;             % set movie time step, in s
%
% Normalize seismograms in some way (after we de-meaned).
% Here we scale seismograms to +- 10. When we use scatter we will plot on a
% scale of +-1, so that we fully saturate the color scale.
%
% How you scale (or don't scale) the seismograms has a big impact on the
% appearance of the visualization. Play aroudn with this.
for i=1:length(dirlist)
    ampmax = max( abs(sacdata(:,i)) ) * 0.1; % use .1 to saturate scale
    sacdata(1:npts(i),i) = sacdata(1:npts(i),i) / ampmax;
end;
fprintf('normalized traces\n');
%
% ADJUST COORDINATES TO MERCATOR PROJECTION (centered on -115 longitude)
% (okay, this is an ugly hack, but it works)
%
lon = lon + 115.0;
lat = asinh(tan(lat*pi/180.))* 180.0 / pi;
shore(:,1) = shore(:,1) + 115.;
shore(:,2) = asinh(tan(shore(:,2)*pi/180.))* 180.0 / pi;
political(:,1) = political(:,1) + 115.;
political(:,2) = asinh(tan(political(:,2)*pi/180.)) * 180.0 / pi;
%
% Okay - ready to start plotting
%
% first, open up a figure and draw the base map
%
fig_h = figure('Position',[1 100 594 714],'Color','white');
hold on
% plot the base map and all parts that don't chnge from frame to frame
% plot(shore(:,1),shore(:,2));
% plot(political(:,1),political(:,2));
% limits for projected coordinates
xlim([-11.0 18.0]);
ylim([30.0 57.0]);
title('Gulf of California, August 3, 2009','FontSize',16);
set(gca,'CLim',[-1.0 1.0]); % set the limits for the call to scatter
set(gca,'Box','on','XTick',[],'YTick',[],'DataAspectRatio',[1 1 1]);
%
% Primary time loop to generate frames of movie
%
% set up a uniform x y grid we will use for contouring - using the
% transformed coordinates
%
 x = -11.0 : 0.5 : 18.0;
 y =  30.0 : 0.5 : 57.0;
[xx yy] = meshgrid (x,y);
colormap('jet');

k = 0;
while tstart < tend
    k = k + 1;
    % extract time time slice from network
    for i=1:length(dirlist)
        tdiff = ( tstart - tbeg(i) ) / dt(i);
        m = int32(tdiff) + 1;  % index into seismogram
        if m<1 || m>npts(i)
            amp(i)= 0.0 ; % seismogram does not overlap target time sample - set value to zero
        else
            amp(i) = sacdata(m,i);
        end;
    end;
    %
    % plot the blinking circles and the time tag
    %
    % generate the mesh of z values that correspond to the amplitudes at
    % this time
    %
    zz = griddata(lon,lat,amp,xx,yy);
    %
    % do the contouring
    [C h1] = contourf(xx,yy,zz);
    set(gca,'CLim',[-10.0 10.0]);
%     h1 = scatter(lon,lat,100,amp,'filled','MarkerEdgeColor','k');
    date_stamp = datestr(tstart/86400,'mm/dd/yyyy HH:MM:SS');
    h2 = text(-11,32,date_stamp,'FontSize',14);
    %
    % plot shore and political on top of the contours, rather than as part
    % of the base map
    %
    h3 = plot(shore(:,1),shore(:,2));
    h4 = plot(political(:,1),political(:,2));

    % Grab the frame, increment the time step, and loop around.
    % First, insert a short pause. I have noticed that sometimes
    % MATLAB grabs the same frame twice if you don't have a short pause.
    drawnow;
    F(k) = getframe(gcf);
    delete(h1);
    delete(h2);
    delete(h3);
    delete(h4);

    tstart = tstart + tstep;
end;
%
% play the movie in matlab (can adjust frame rate if desired - see help file)
% movie(F);
%
% write the movie to an avi file (can adjust frame rate if desired - see
% help file)
% movie2avi(F,'movie_name');
