Contents of file named script.m

% read a set of traces into an array of SAC structures. 
RS = msac_mread('*.BHE') ;

% build normalised amplitude, time and distance matrices.
for itr = 1:length(RS)
   A(:,itr) = RS(itr).x1./(RS(itr).depmax-RS(itr).depmin) ;
   D(:,itr) = ones(1,RS(itr).npts).*RS(itr).gcarc ;
   T(:,itr) = (0:RS(itr).npts-1).*RS(itr).delta + RS(itr).b ;
end

% sort matrices into order of increasing epicentral distance.
[~,ind] = sort([RS(:).gcarc]) ;
AS = A(:,ind) ; TS = T(:,ind) ; DS = D(:,ind) ;

% plot as greyscale surface plot.
surf(DS,TS,AS,'LineStyle','none')
colormap(gray) ; caxis([-0.75 0.75]); axis tight
xlabel('Epicentral distance (deg)') ; ylabel('Time (s)')
shading interp ; view(0,90) ;