% 
% Simple MATLAB script to plot global stacks available here:
% http://www.iris.edu/ds/products/globalstacks_supplement
% Alex Hutko, IRIS DMC
% 2014.09.02
%
% This example uses:
% http://www.iris.edu/files/products/globalstacks/GlobalStacksFiles_BHZ_1_sec.zip
%

clear all;

TITLE = 'Example stack figure'
FIGURENAME = 'ExampleStackFigure'

%------------------------ read in files

PrefixString = 'GlobalStacks.BHZ.1.sec.dist'
for ibin=1:360;
   ibin
   infile=sprintf('%s.%d',PrefixString,ibin);
   a=load(infile);
   dist(ibin)=ibin;
   time=a(:,1);
   STA=a(1,2);
   NTraces(ibin)=a(1,3);
   NRejected(ibin)=a(1,4);
   StackRaw=a(:,5);
   StackRawN2=a(:,6);
   StackSTALTA=a(:,7);
   StackSTALTAN2=a(:,8);
   StackEnvelope=a(:,9);

   ThingToPlot= sqrt(abs(StackSTALTAN2 - mean(StackSTALTAN2)));   %--- this is where you get creative and clean up the image

   z(1:length(time),ibin)=ThingToPlot;
end

%------------------------ set up x and y axes

x=dist;
y=time./60;

%------------------------ plot main figure


fig=figure('units','pixels','outerposition',[20 20 500 600 ],'paperpositionmode','auto','renderer','zbuffer');    %   left up width height
clf
subplot('Position',[ 0.13 0.21 0.82 0.71 ]);  %   left up width height
set( gca,'Fontsize',12,'FontName','Times')

pcolor(x,y,z)   %--- magic happens here

shading flat
titlename = sprintf(' %s ',TITLE)
title(titlename)
xlabel('Distance (deg)')
ylabel('Time (minutes)')
set(gca,'Xticklabel',[ ])
colorbar off
axis([ 0 360 0 90 ])

%------------------------ play with colors

boneinv(length(colormap(bone)):-1:1,1:3)=colormap(bone);
colormap(boneinv)
%colormap(jet)
%colormap(copper)

SaturationFactor=0.3;

caxis([ 0 max(max(z))*SaturationFactor ])
ylabel(' Time (minutes) ')

%------------------------ plot histogram of traces

subplot('Position',[ 0.13 0.09 0.82 0.11 ]);  %   left up width height
set( gca,'Fontsize',12,'FontName','Times')

bar(0:0.5:max(dist)/2-0.5,NTraces,'facecolor',[0.5 0.6 0.6])

ylabel(' # of traces ')
xlabel(' Distance (degrees) ')
axis([ 0 180 0 10000 ])
set( gca,'Fontsize',12,'FontName','Times')
sumcount = sum(NTraces);
txttotalcount=sprintf(' %4d traces ',sumcount)
text(130,7000,txttotalcount,'backgroundcolor','w','Fontsize',12,'FontName','Times')

%------------------------ print to .png file

printname = sprintf('%s.png',FIGURENAME)
print('-dpng','-r250',printname)

