%Script for generating the 'Stack Traces' from the original
%  array data in Seismic Handler format.
%
% Requires:      rqhd.m  - read seismic handler header files.
%                rqbn.m  - read seismic handler binary files.
%
%           filelist     - file containing a list of seismic
%                          handler binary files to process.
%          
%           P_ttimes.xy  - file containing approximate first
%                          arrival traveltimes.
%                           

clear all

% enter input variables
%-------------------------------------------------------------------------
nevents     =324;   %Total number of events in temp_binlist


%Code requires P_ttimes.xy
%-------------------------------------------------------------------------
load('P_ttimes.xy');
fidd=fopen(['filelist']);             %list of data files to process

while 1
infile=fgetl(fidd)                    %loop through all files
if ~ischar(infile), break, end

tic

%read in data
%-------------------------------------------------------------------------
data=rqbn(infile);
tlength=length(data(1,:));
delta=rqhd(infile,'R000',1);
tt=size(data);
for k=1:tt(1,1)
  slow(k)=str2num(rqhd(infile,'S000:S',k));
end

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


%get event times  
%-------------------------------------------------------------------------
%origin time
ot1=rqhd(infile,'S024',1);
hour=str2num(ot1(13:14))*3600;
minutes=str2num(ot1(16:17))*60;
sec=str2num(ot1(19:length(ot1)));
ot=hour+minutes+sec;

%begin time
bt1=rqhd(infile,'S021',1);
hour=str2num(bt1(13:14))*3600;
minutes=str2num(bt1(16:17))*60;
sec=str2num(bt1(19:length(ot1)));
bt=hour+minutes+sec;

%start time of trace in seconds
at=bt-ot;


%determine time values
%-------------------------------------------------------------------------
xvals=linspace(at,at+tlength*delta,tlength);


%collapse original data to stack trace
%-------------------------------------------------------------------------

[dmin, mindex] = min(data);
[dmax, maxdex] = max(data);

% saved variable is orx 3-columns (time,amplitude,slowness) 
%
for i=1:tlength
  if dmax(i) > abs(dmin(i))
    orx(i,2) = dmax(i);
    orx(i,3) = slow(maxdex(i));
  else
    orx(i,2) = dmin(i);
    orx(i,3) = slow(mindex(i));
  end
end
orx(:,1)=xvals';

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

%zero out values in 'orx' that are below the average noise level
IB = find(orx(:,2) < Alevel);
for kk=1:length(IB)
  orx(IB(kk),2) = 0.0;
  orx(IB(kk),3) = 0.0;
end

figure(1)
plot(orx(:,1),orx(:,2))
title('matlab variable orx')
xlabel('time (sec)')
ylabel('non-normalized amplitude')
grid on

%save variables in mat file
%-------------------------------------------------------------------------
outname=strcat(infile,'_prx');
evdp=rqhd(infile,'R014',1);  %event depth
lon=rqhd(infile,'R017',1);   %event longitude
lat=rqhd(infile,'R016',1);   %event latitude
azi=rqhd(infile,'R012',1);   %azimuth
gcarc=rqhd(infile,'R011',1); %distance

save(outname,'orx','evdp','lon','lat','azi','gcarc') 

toc

clear outname evdp lon lat azi gcarc
clear orx data tlength delta slow ot1 hour minutes sec ot bt1  
clear bt at xvals dmin mindex dmax maxdex i infile 


end
fclose(fidd);