%RAW_AMPLITUDE_PLOT
%
%Produce Amplitude Plot and write amplitudes/slownesses to binary
%  x,y,v files readable with GMT
%
%  Input:  needs a file 'matfile_list' which lists the files
%          to use in the stacking procedure.
%
%  Dependencies: none
%
%  By, Michael Thorne (mthorne@gi.alaska.edu)
%
%  Last Update:  Sept. 2005

clear all

%Set options
%**********************************************************************
nevents  =1796;  %number of events in 'matfile_list'
ddist    =0.5;   %distance bin size (degrees)
dt       =0.5;   %time bin size (sec)
mindist  =4;     %minimum distance of plot (deg)
maxdist  =180;   %maximum distance of plot(deg)
mintime  =0;     %minimum time of plot (sec) 
maxtime  =2000;  %maximum time of plot (sec)
nonlin   =0.5;   %nonlinear scale (~0-1.0; 1=no scaling)
noise    =1;     %0=don't whiteout noise, 1=whiteout noise
nscale   =1;     %set to zero all amplitudes <= nscale*(average noise)
fil      =0;     %0=no filtering, 1=filter envelopes
lowp     =0.5;   %low pass filter cutoff frequency [Hz]
poles    =2;     %number of low pass filter poles
passes   =1;     %number of low pass filter passes
writeout =0;     %write binary files (0=no, 1=yes)
%**********************************************************************


load('P_ttimes.xy');  %Approximate P-wave arrival times

%Set up grid
%**********************************************************************
tgrid=0:dt:4000;
dgrid=0:ddist:180;
dsize=size(dgrid);
tsize=size(tgrid);
[dd,tt]=meshgrid(dgrid,tgrid);
%**********************************************************************

%Initialize arrays
%**********************************************************************
msize=size(dd);
super(1:msize(1,1),1:msize(1,2))= 0.0;       %amplitude array
navg(1:msize(1,1),1:msize(1,2))= 1;          %no. traces per distance
slow_avg(1:msize(1,1),1:msize(1,2))= 0.0;    %average slowness
%**********************************************************************

%Begin loop through all files (distance bin loop)
%**********************************************************************
fid=fopen(['matfile_list']);
for event=1:nevents
 textmes=[num2str(event),'/',num2str(nevents)];
 disp(textmes)                                  %display program status

  infile = deblank(fgets(fid));
  load(infile,'-mat');                          %load matfile

  distance = gcarc;
  dbin=round(distance/ddist);                   %assign distance bin

  %Whiteout below noise level
  %********************************************************************
  if (noise == 1)

    %determine approximate P-arrival time
    %********************************************************************
    for i=1:length(P_ttimes(:,1))
      if (round(distance) == P_ttimes(i,2))
        ptime = P_ttimes(i,1);
      end
    end

    %determine average noise level before direct P-arrival
    IA = find(orx(:,1) < (ptime-20));
    TA = sum(orx(IA,2));
    Alevel = TA/length(IA);

    %zero out values in 'orx' that are below the average noise level
    IB = find(orx(:,2) <= nscale*Alevel);
    orx(IB,2) = 0.0;
    orx(IB,3) = 0.0;

    clear IB TA Alevel IA
  end
  %********************************************************************
  
  %Grab time,amplitude,slowness info from matfile
  %********************************************************************
  sm_size = size(orx);
  xvals=orx(:,1);                               %time axis
  scale=max(abs(orx(:,2)));                     %maxval amplitudes
  yvals=orx(:,2)/scale;                         %normalized amplitudes
  zvals=orx(:,3);                               %slowness

  %Low Pass filter trace
  %********************************************************************
  if (fil == 1)
    %1/2 Sample rate [Hz]
    delta = (1.0/(xvals(2) - xvals(1)))/2.0;    
    %Butterworth lowpass filter
    %[ButterB,ButterA]= butter(poles,lowp/delta,'low'); 
    [ButterB,ButterA]=cheby2(poles,20,lowp/delta,'low');
    %tfilt = filter(ButterB,ButterA,yvals);
    %yvals = tfilt;
    for kkk=1:passes
      yvals = filter(ButterB,ButterA,yvals);
    end
  end

  %loop through time bins
  %********************************************************************
  for j=1:sm_size(1,1)
    time=xvals(j);
    addval=yvals(j);
    baddval=zvals(j);
    tbin=round(time/dt); 
    if (tbin > 0)
      super(tbin,dbin) = super(tbin,dbin) + addval;
      slow_avg(tbin,dbin) = slow_avg(tbin,dbin) + baddval;
      navg(tbin,dbin)=navg(tbin,dbin)+1;
    end
  end

  %Clear matlab variable before reading in next event
  %********************************************************************
  clear dif ff tbin time addval
  clear gcarc azi distance evdp lat lon orx 

end

%readjust navg
%**********************************************************************
for dbin=1:length(dgrid)
  for tbin=1:length(tgrid)
    if (navg(tbin,dbin) > 1)
      navg(tbin,dbin) = navg(tbin,dbin) - 1;
    end
  end
end

%average the values within each bin
%**********************************************************************
super=super./navg;

%scale entire range of amplitudes by maximum 
%**********************************************************************
scale=max(max(super));
super=super./scale;

%non-lin scaling
%**********************************************************************
super=(abs(super)).^nonlin;

%average slowness
%**********************************************************************
slow_avg=slow_avg./navg;

%set unsampled regions in slowness to low value
%**********************************************************************
for dbin=1:length(dgrid)
  for tbin=1:length(tgrid)
    if (slow_avg(tbin,dbin) == 0)
      slow_avg(tbin,dbin) = -999.0;
    end
  end
end

if writeout == 1

%write out amplitudes to file readable with GMT
%**********************************************************************
fid = fopen('superfile.bin','w','ieee-le')
for dbin=1:length(dgrid)
  for tbin=1:length(tgrid)
    outp(1) = dgrid(dbin);
    outp(2) = tgrid(tbin);
    outp(3) = super(tbin,dbin);
    fwrite(fid,outp,'single');
  end
end
fclose(fid)

%write out slowness to file readable with GMT
%**********************************************************************
fid = fopen('slowness.bin','w','ieee-le')
for dbin=1:length(dgrid)
  for tbin=1:length(tgrid)
    outp(1) = dgrid(dbin);
    outp(2) = tgrid(tbin);
    outp(3) = slow_avg(tbin,dbin);
    fwrite(fid,outp,'single');
  end
end
fclose(fid)

end

%Plot Amplitudes
%**********************************************************************
figure(1)
colormap(hot)
set(gcf,'PaperOrientation','landscape','PaperPosition', ...
    [.25 .25 10.5 8],'PaperType','usletter');
imagesc(dgrid,tgrid,super)
set(gca,'YDir','normal')
caxis([0 1])
axis([mindist maxdist mintime maxtime])
xlabel('Distance (deg)')
ylabel('Time (sec)')
title('Amplitude')
colorbar
%**********************************************************************