%Script for generating the 'Stack Traces' from the original
%  array data in raw format.
%
% Requires: filelist - file containing a list of files
%                          to process 
%          
%           P_ttimes.xy  - file containing approximate first
%                          arrival traveltimes.
%                           

clear all


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

while 1                         %loop through all files in 'filelist'
infile=fgetl(fidd)              %    returns -1 
if ~ischar(infile), break, end

tic

%read in data
%-------------------------------------------------------------------------
fid2=fopen(infile);

%read header info
delta = str2num(fgetl(fid2));    %dt    - sampling - time increment (sec)
gcarc = str2num(fgetl(fid2));    %gcarc - great circle distance (deg)
ot1   = str2num(fgetl(fid2));    %ot1   - event origin time (sec)
bt1   = str2num(fgetl(fid2));    %bt1   - trace begin time (sec)
evdp  = str2num(fgetl(fid2));    %evdp  - event depth (km)
lon   = str2num(fgetl(fid2));    %lon   - event longitude (deg)
lat   = str2num(fgetl(fid2));    %lat   - event latitude (deg)
azi   = str2num(fgetl(fid2));    %azi   - event/rec azimuth (deg)

%read slownesses and data
counter = 1;
while 1
  tline = fgetl(fid2);
  if ~ischar(tline), break, end
  temp = str2num(tline);
  slow(counter) = temp(1,1);
  iend = size(temp);
  data(counter,:) = temp(1,2:iend(2));
  counter = counter + 1;
end
fclose(fid2);
clear temp

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

%get event times  
%-------------------------------------------------------------------------
%start time of trace in seconds
at=bt1-ot1;


%determine time values
%-------------------------------------------------------------------------
tlength = length(data(1,:));
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

%pause

%save variables in mat file
%-------------------------------------------------------------------------
outname=strcat(infile,'_prx');

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

toc

clear outname evdp lon lat azi gcarc orx data tlength delta slow ot1 bt1  
clear at xvals dmin mindex dmax maxdex i infile IA TA Alevel IB

end
fclose(fidd);