% This script is a modified version of the first MATLAB script, which factors the HIPR weights into
% the data stacks.
%
% The code expects the input list of events to have three additional columns (in addition to the event
% name). The second column lists either "B" or "O", indicating which list (best or other) the event came 
% from.  The third column lists either "1" (meaning to keep/run event) or "0" (meaning to skip the event).  
% The last column indicates the position/index of the event in the corresponding event-weight input file.  
% Done to ensure the correct weight values are matched up with the correct event.
%
% Samantha Hansen, 02/20/20
%

function varargout = PcP_GUI(varargin)
%  PcP_GUI MATLAB code for PcP_GUI.fig
%  Created by Samantha Hansen
%  Last modified: June 20, 2020

%% BEGIN INITIALIZATION - DO NOT EDIT THIS SECTION %%
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @PcP_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @PcP_GUI_OutputFcn, ...
                   'gui_LayoutFcn',  [], ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
  gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
  [varargout{1:output}] = gui_mainfcn(gui_State, varargin{:});
else
  gui_mainfcn(gui_State, varargin{:});
end

%% Function to setup PcP_GUI and make it visible
function PcP_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
  % This function has no output args, see OutputFcn.

  % Choose default command line output for PcP_GUI
  handles.output = hObject;

  % Update handles structure
  guidata(hObject, handles);

%% Function outputs are returned to the command line
function varargout = PcP_GUI_OutputFcn(hObject, eventdata, handles)

  % Get default command line output from handles structure
  varargout{1} = handles.output;

%% END OF INITIALIZATION - DO NOT EDIT THE ABOVE SECTION%%


%% START FUNCTION
%% Function executes when "Start" button is pressed

function Start_Callback(hObject, eventdata, handles) 
  
  % Read in data of interest from user-specified list
  readme = importdata('ISlist');
  handles.list = readme.textdata(:,1);
  handles.BorO = readme.textdata(:,2);
  handles.KeepDump = readme.data(:,1);
  handles.WeightIndex = readme.data(:,2);
    
  % Read in station and event weight files; add station names in first column of station weight matrix
  % (note: event weights in same order so don't need to do for those matricies)
  load './stationWeights'; 
  load './eventWeights_P'; 
  load './eventWeights_PcP'; 
  load './eventWeights_P_otherEvts'; 
  load './eventWeights_PcP_otherEvts'; 
  stations(1,1) = {'APRL'};
  stations(2,1) = {'BEBP'};
  stations(3,1) = {'CASY'};
  stations(4,1) = {'DONT'};
  stations(5,1) = {'DUBY'};
  stations(6,1) = {'FOOT'};
  stations(7,1) = {'GRAW'};
  stations(8,1) = {'KNYN'};
  stations(9,1) = {'LEON'};
  stations(10,1) = {'MICH'};
  stations(11,1) = {'RAPH'};
  stations(12,1) = {'RKST'};
  stations(13,1) = {'SAMH'};
  stations(14,1) = {'SHRD'};
  stations(15,1) = {'SPLN'};
  for i=1:length(stationWeights)
    SW(i,1) = stations(i,1);
    SW(i,2) = {stationWeights(i,1)};
  end

  % For each event in the list, the process starts here
  SkipEvt = 0;
  for counter = 1:length(handles.list)

    handles.k = counter;

    if (handles.k > 1)
      % Clear fields for next iteration
      handles = rmfield(handles,'eventname');
      handles = rmfield(handles,'eventfolder');
      handles = rmfield(handles,'path_to_P_data');
      handles = rmfield(handles,'stationW');
      handles = rmfield(handles,'evtWP');
      handles = rmfield(handles,'evtWPcP');
      set(handles.text1, 'String','');
    end
    if (handles.k > 1 && SkipEvt == 0)
      handles = rmfield(handles,'dataP');
      handles = rmfield(handles,'dataPcPorig');
      handles = rmfield(handles,'dataPcut');
      handles = rmfield(handles,'dataPcPcut');
      handles = rmfield(handles,'sumShiftP'); 
      handles = rmfield(handles,'sumShiftPcP'); 
      handles = rmfield(handles,'oldPstack');
      handles = rmfield(handles,'dataPcopy');
      handles = rmfield(handles,'timeShiftP');
      handles = rmfield(handles,'coefP');
      handles = rmfield(handles,'weightPstack');
      handles = rmfield(handles,'snrP');
      handles = rmfield(handles,'dataPcP');
      handles = rmfield(handles,'oldPcPstack');
      handles = rmfield(handles,'dataPcPcopy');
      handles = rmfield(handles,'timeShiftPcP');
      handles = rmfield(handles,'coefPcP');
      handles = rmfield(handles,'weightPcPstack');
      handles = rmfield(handles,'snrPcP');
      handles = rmfield(handles,'weightPcPstackPad');
      set(handles.text2, 'String','');
      set(handles.text3, 'String','');
      set(handles.text17, 'String','');
      set(handles.text20, 'String','');
      cla(handles.axes1,'reset');
      cla(handles.axes2,'reset');
      cla(handles.axes3,'reset');
      cla(handles.axes6,'reset');
      cla(handles.axes7,'reset');
      cla(handles.axes8,'reset');
      cla(handles.axes12,'reset');
    end
    % Reset SkipEvt flag if needed
    SkipEvt = 0;

    % Read in specified event from list 
    handles.eventname = handles.list(handles.k);
    set(handles.text1, 'String', handles.eventname);

    % Specify path to data
    handles.eventfolder = strcat('./PcP_Data/',handles.list(handles.k),'/FILT_0.5-1.5/');
    handles.path_to_P_data = strcat('./PcP_Data/',handles.list(handles.k),'/FILT_0.5-1.5/*.Pf');

    % Specify cut limits, signal window sizes, correlation window size, and max time shift allowance
    cutMin = -25;		% Minimum time limit for cut records 
    cutMax = 25;		% Maximum time limit for cut records 
    signalStart = -3;		% Number of secs prior to reference time on individual records to start signal window
    signalEnd = 10;		% Number of secs after the reference time on individual records to end signal window 
    noiseStart = -25;		% Number of secs prior to reference time on individual records to start noise window 
    noiseEnd = -10;		% Number of secs prior to reference time on individual records to end noise window
    signalStartStack = -3;	% Number of secs prior to max stack amplitude to start signal window 
    signalEndStack = 7;		% Number of secs after the max stack amplitude to end signal window 
    noiseStartStack = -20;	% Number of secs prior to max stack amplitude to start noise window 
    noiseEndStack = -10;	% Number of secs prior to max stack amplitude to end noise window 
    limitTimeMax = 10;	        % Number of secs +/- around predicted time to search for P or PcP arrival (upper limit)
    limitTimeMin = 5;		% Number of secs +/- around predicted time to search for P or PcP arrival (lower limit)
    windowSizeT = 5;		% +/- time to isolated signal for cross-correlation
    maxShift = 4;		% max time shift allowance 
    codaAllow = 15;		% Number of seconds "buffer" to allow for P-wave delays and coda (when checking position vs. PcP noise window)
    sumPhaseTimeDiff = 0;	% Initialize value to find average PcP-P time (from PREM)

    % Read in P and PcP data, sort by distance, and cut P records
    clearvars AlldataP distance sortDist WP WPcP
    AlldataP = readseizmo(handles.path_to_P_data);
    for i=1:length(AlldataP)
       distance(i,1:2)={getheader(AlldataP(i),'DIST') getheader(AlldataP(i),'KSTNM')};
    end
    sortDist = sortrows(distance,1);

    % Get associated event weights from correct position in weighting matrices
    if (strcmp(handles.BorO(handles.k),'B'))
      WP(:,1) = eventWeights_P(:,handles.WeightIndex(handles.k));
      WPcP(:,1) = eventWeights_PcP(:,handles.WeightIndex(handles.k));
    else
      WP(:,1) = eventWeights_P_otherEvts(:,handles.WeightIndex(handles.k));
      WPcP(:,1) = eventWeights_PcP_otherEvts(:,handles.WeightIndex(handles.k));
    end

    % Reorder weights to match this event (sorted by distance)
    clearvars Cell staNameE Cell2 staNameM 
    for i=1:length(sortDist)
      Cell = sortDist(i,2);
      staNameE = Cell{1};
      for j=1:length(SW)
        Cell2 = SW(j,1);
        staNameM = Cell2{1};
        if(strcmp(staNameM,staNameE))
          % Found a matching pair; reorder weights to match distance sorting
          handles.stationW(i,1) = cell2mat(SW(j,2));
          handles.evtWP(i,1) = WP(j,1);
          handles.evtWPcP(i,1) = WPcP(j,1);
        end
      end
    end

    % Put in check here; if all event weights for P and/or PcP are zero OR if the "KeepDump" flag
    % is set to zero, then skip event
    clear sumWP sumWPcP
    sumWP = 0;
    sumWPcP = 0;
    for i=1:length(handles.evtWP)
      sumWP = sumWP + handles.evtWP(i,1);
      sumWPcP = sumWPcP + handles.evtWPcP(i,1);
    end
    if (sumWP == 0 || sumWPcP == 0 || handles.KeepDump(handles.k) == 0)
      % Have a case where all P and/or PcP event weights are zero; want to skip this event
      SkipEvt = 1;
      continue
    end

    clearvars Cell staName PstaName PcPstaName maxampP sampRate pickTimes timeP PREMavg timePcP SNRweightP SNRweightPcP
    for i=1:length(sortDist)
      Cell = sortDist(i,2);
      staName = Cell{1};
      PstaName = strcat(handles.eventfolder,staName,'.Pf');
      PcPstaName = strcat(handles.eventfolder,staName,'.P.PcPf');
     
      % Full P and PcP records (not cut); reading them in in this order sorts them by distance
      handles.dataP(i) = readseizmo(PstaName);
      handles.dataPcPorig(i) = readseizmo(PcPstaName);
    
      % Save original (PREM-predicted) pick times for other phases (T6 = pP, T7 = sP) relative to both P and PcP and adjusted so that the "main" phase (P or PcP) is at zero
      pickTimes(i,1) = ((getheader(handles.dataP(i),'T6'))-cutMax) - ((getheader(handles.dataP(i),'T1')) - cutMax);
      pickTimes(i,2) = ((getheader(handles.dataP(i),'T7'))-cutMax) - ((getheader(handles.dataP(i),'T1')) - cutMax);
      pickTimes(i,3) = ((getheader(handles.dataPcPorig(i),'T6'))-cutMax) - ((getheader(handles.dataPcPorig(i),'T4')) - cutMax);
      pickTimes(i,4) = ((getheader(handles.dataPcPorig(i),'T7'))-cutMax) - ((getheader(handles.dataPcPorig(i),'T4')) - cutMax);

      % Add up PcP-P time (based on PREM) so can find average
      sumPhaseTimeDiff = sumPhaseTimeDiff + ((getheader(handles.dataPcPorig(i),'T4')) - (getheader(handles.dataP(i),'T1')));

      % Cut P records using limits specified above and normalize amplitudes 
      handles.dataPcut(i) = cut(handles.dataP(i),'t1',cutMin,'t1',cutMax);
      maxampP = max(handles.dataPcut(i).dep);
      handles.dataPcut(i).dep = (handles.dataPcut(i).dep) / maxampP;
      sampRate = 1 / getheader(handles.dataP(i),'DELTA');

      % Define time axis (used for all P records); adjust so PREM P time is at zero
      timeP = ((1:length(handles.dataPcut(i).dep))/sampRate) - cutMax;
     
      % Do similar thing for PcP but won't normalize amplitudes right now (just setting up time); again, adjusted so PREM PcP time is at zero
      handles.dataPcPcut(i) = cut(handles.dataPcPorig(i),'t4',cutMin,'t4',cutMax);
      timePcP = ((1:length(handles.dataPcPcut(i).dep))/sampRate) - cutMax;
      handles = rmfield(handles,'dataPcPcut');
 
      % Initialize  matricies to keep track of time shifts (both total and absolute) and initialize matrices for the SNR weights
      handles.sumShiftP(i) = 0;
      handles.sumShiftPcP(i) = 0;
      SNRweightP(i) = 1;
      SNRweightPcP(i) = 1;
    end
      
    % Convert windowSizeT (in time) to points, using the data sampling rate
    windowSizeP = int32(windowSizeT * sampRate);

    % Compute and report average PcP-P time (based on PREM)
    PREMavg = sumPhaseTimeDiff / (length(sortDist));
    set(handles.text17, 'String',(num2str(PREMavg)));
 
    % Display the event info at the top of the GUI window
    depth = getheader(handles.dataP(1),'EVDP');
    depth = num2str(depth);
    lat = getheader(handles.dataP(1),'EVLA');
    lat = num2str(lat);
    lon = getheader(handles.dataP(1),'EVLO');
    lon = num2str(lon);
    mag = getheader(handles.dataP(1),'MAG');
    mag = num2str(mag);
    textLabel = strcat('Depth:',32,depth,32,'Lat:',32,lat,32,'Lon:',32,lon,32,'Mag:',32,mag);
    set(handles.text2, 'String',textLabel);
   
    % First make linear P-wave stack given predicted P alignment (no weighting); normalize the stack so max amplitude = one
    clearvars pStack0 count MT limitTimeP
    pStack0 = 0;
    count = 0;
    for i = 1:length(handles.dataPcut)
       pStack0 = pStack0 + handles.dataPcut(i).dep;
       count = count + 1;
    end
    handles.oldPstack = pStack0 / count;
    handles.oldPstack = (handles.oldPstack) / (max(handles.oldPstack));
   
    % Find the time of the max amplitude in relation to the PREM predicted P time
    % use that time different to control search window
    MT = ((find(handles.oldPstack == max(handles.oldPstack))) / sampRate) - cutMax; 
    if (MT > limitTimeMax)
      limitTimeP = int32(limitTimeMin * sampRate);
    else
      limitTimeP = int32(limitTimeMax * sampRate);
    end 
    
    % In some cases, just defining P as the largest amplitude on the stack doesn't work well - need to give the code
    % some timing guidance.  Assume that P arrives within a +/- limitTime of the predicted arrival time (zero on the
    % plots).  Find the max amplitude within that window; extract a portion of the P stack +/- windowSizeT around that
    % max amplitude; taper and pad with zeros so length is the same as stack.  Will use similar approach with PcP below.
    clearvars searchWinStart searchWinEnd maxampP posA startWin endWin extractPsignal taperedP ZF ZB TP
    searchWinStart = ((length(handles.oldPstack))/2)-limitTimeP;
    searchWinEnd = ((length(handles.oldPstack))/2)+limitTimeP;
    maxampP = 0;
    for (i = 1:length(handles.oldPstack))
      if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.oldPstack(i,1) > maxampP))
        maxampP = handles.oldPstack(i,1);
        posA = i;
      end
    end
    startWin = posA - windowSizeP;
    endWin = posA + windowSizeP;
    extractPstack = handles.oldPstack(startWin:endWin);
    taperedPstack = (extractPstack .* tukeywin((length(extractPstack)),0.25));
    ZF = (zeros([1,(startWin-1)]))';
    ZB = (zeros([1,(length(handles.oldPstack)-endWin)]))';
    TPstack = [ZF;taperedPstack;ZB];

    % Cut and taper the individual P records the same as the stack (use same startWin and endWin)
    clearvars extractPsignal taperPsignal ZF ZB TPsignal acor lag I acorr t1 newT1 maxampP maxCheck saveMaxAmp
    handles.dataPcopy = handles.dataP;
    for i = 1:length(handles.dataPcut)
      extractPsignal = handles.dataPcut(i).dep(startWin:endWin);
      taperPsignal = (extractPsignal .* tukeywin((length(extractPsignal)),0.25));
      ZF = (zeros([1,(startWin-1)]))';
      ZB = (zeros([1,(length(handles.dataPcut(i).dep)-endWin)]))';
      TPsignal = [ZF;taperPsignal;ZB];

      % Cross-correlate the extracted/tapered P stack with each extracted/tapered individual P record
      % to determine a time shift and the associated CCC values
      [acor,lag] = xcorr(TPsignal,TPstack);
      [~,I] = max(acor);
      lagDiff(i) = lag(I);
      handles.timeShiftP(i) = lagDiff(i) / sampRate;
      acorr = xcorr(TPsignal,TPstack,'coeff');
      handles.coefP(i) = max(abs(acorr));
      handles.sumShiftP(i) = handles.sumShiftP(i) + handles.timeShiftP(i);
    
      % Apply the correlated time shift to adjust the t1 marker, capping the shift at +/- maxShift
      if ((abs(handles.sumShiftP(i))) >= maxShift)
        if (handles.sumShiftP(i) < 0)
          handles.sumShiftP(i) = (-1 * maxShift);
        else
          handles.sumShiftP(i) = maxShift;
        end
      end

      t1(i) = getheader(handles.dataPcopy(i),'t1');
      newT1(i) = t1(i) + handles.sumShiftP(i);
      handles.dataPcopy(i) = changeheader(handles.dataPcopy(i),'t1',newT1(i));
      handles.dataPcut(i) = cut(handles.dataPcopy(i),'t1',cutMin,'t1',cutMax);

      % Before normalize waveform amps, save their max amplitude
      maxCheck = 0;
      for (j = 1:length(handles.dataPcut(i).dep))
        if ((j >= searchWinStart) && (j <= searchWinEnd) && (handles.dataPcut(i).dep(j) > maxCheck))
          maxCheck = handles.dataPcut(i).dep(j); 
        end
      end
      saveMaxAmpP(i) = maxCheck;
 
      % Normalize waveform amps before moving forward
      maxampP = max(handles.dataPcut(i).dep);
      handles.dataPcut(i).dep = (handles.dataPcut(i).dep) / maxampP;
    end

    % At this point, all data has full weight.  Setup loop to determine SNRs, which will be combined with CCC values
    % to weight each waveform.  The weighted waveforms are combined into a weighted (normalized) stack.  The loop will 
    % continue as long as the new stacks are getting better and until they are only minimally different from the prior
    % one.  If during this process a record hits the maximum time shift allowance, it is given an SNR weight of zero 
    % to remove it from stacks and later analyses.
    clearvars compare compareOld
    compare = 0.01;
    compareOld = 0;

    while (compare <= 0.995) 
      if (compareOld < compare)
        clearvars count 
        handles.weightPstack = 0;
        count = 0;

        for i = 1:length(handles.dataPcut)

           clearvars signalPsamp noisePsamp signalP noiseP acor weightP

           % Compute the signal-to-noise ratio for each individual waveform
           signalPsamp = cut(handles.dataPcut(i),'t1',signalStart,'t1',signalEnd);
           noisePsamp = cut(handles.dataPcut(i),'t1',noiseStart,'t1',noiseEnd);
           signalP = mean(abs(signalPsamp.dep));
           noiseP = mean(abs(noisePsamp.dep));
           handles.snrP(i) = signalP / noiseP;
    
           % Setup ramp function for SNR weighting (between specified max and min limits)
           if ((handles.snrP(i) >= 10) && (SNRweightP(i) ~= 0))
             SNRweightP(i) = 1;
           elseif ((handles.snrP(i) < 10) && (handles.snrP(i) > 1) && (SNRweightP(i) ~= 0))
	     % Linear function between two extreme values; NOTE: if change the weighting limits
             % will need to change this equation too 
             SNRweightP(i) = (0.1056*handles.snrP(i)) - 0.0556;       
           elseif ((handles.snrP(i) <= 1) && (SNRweightP(i) ~= 0))
             SNRweightP(i) = 0.05;
           end
 
           weightP = handles.dataPcut(i).dep * SNRweightP(i) * (handles.coefP(i)^2) * handles.stationW(i,1) * handles.evtWP(i,1);
           if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0) 
             handles.weightPstack = handles.weightPstack + weightP;
             count = count + 1;     
           end
        end
        
        % Compute weighted P stack, normalizing its max amplitude to one
        handles.weightPstack = handles.weightPstack / count;    
        handles.weightPstack = (handles.weightPstack) / (max(handles.weightPstack));

        % Compare the newest weighted stack to the previous stack
        acor = xcorr(handles.oldPstack,handles.weightPstack,'coeff');
        compareOld = compare;
        compare = max(acor);

        % Check similarity between stacks to see if need to keep going
        if (compare <= 0.995)
   
          % Like above, extract P signal from stack, taper, and use this for the next cross-correlation 
          clearvars searchWinStart searchWinEnd maxampP posA startWin endWin extractPstack taperedPstack ZF ZB TPstack 
          searchWinStart = ((length(handles.weightPstack))/2)-limitTimeP;
          searchWinEnd = ((length(handles.weightPstack))/2)+limitTimeP;
          maxampP = 0;
          for (i = 1:length(handles.weightPstack))
            if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPstack(i,1) > maxampP))
              maxampP = handles.weightPstack(i,1);
              posA = i;
            end
          end
          startWin = posA - windowSizeP;
          endWin = posA + windowSizeP;
          extractPstack = handles.weightPstack(startWin:endWin);
          taperedPstack = (extractPstack .* tukeywin((length(extractPstack)),0.25));
          ZF = (zeros([1,(startWin-1)]))';
          ZB = (zeros([1,(length(handles.weightPstack)-endWin)]))';
          TPstack = [ZF;taperedPstack;ZB];
    
          % Cut and taper the individual P records the same as the stack (use same startWin and endWin)
          clearvars extractPsignal taperPsignal ZF ZB TPsignal acor lag I acorr t1 newT1 maxampP maxCheck saveMaxAmpP
          handles = rmfield(handles,'dataPcopy');
          handles.dataPcopy = handles.dataP;
          for i = 1:length(handles.dataPcut)
            extractPsignal = handles.dataPcut(i).dep(startWin:endWin);
            taperPsignal = (extractPsignal .* tukeywin((length(extractPsignal)),0.25));
            ZF = (zeros([1,(startWin-1)]))';
            ZB = (zeros([1,(length(handles.dataPcut(i).dep)-endWin)]))';
            TPsignal = [ZF;taperPsignal;ZB];

            % Cross-correlate the extracted/tapered P stack with each extracted/tapered individual P record
            % to determine a time shift and the associated CCC values
            [acor,lag] = xcorr(TPsignal,TPstack);
            [~,I] = max(acor);
            lagDiff(i) = lag(I);
            handles.timeShiftP(i) = lagDiff(i) / sampRate;
            acorr = xcorr(TPsignal,TPstack,'coeff');
            handles.coefP(i) = max(abs(acorr));
            
            % Add this new time shift to the shift sum (both total and absolute)
            handles.sumShiftP(i) = handles.sumShiftP(i) + handles.timeShiftP(i);

            % Like above, time shift is limited to +/- maxShift
            if ((abs(handles.sumShiftP(i))) >= maxShift)
              if (handles.sumShiftP(i) < 0)
                handles.sumShiftP(i) = (-1 * maxShift);
              else
                handles.sumShiftP(i) = maxShift;
              end
              % If the time shift hits the maximum shift allowance, set SNR weight = 0
              SNRweightP(i) = 0;
            else
              if (SNRweightP(i) == 0)
                 % Was assigned a weight of zero on last iteration because hit max time shift but it has
                 % move away from that limit now.  So reassign its weight using the same scheme as above
                 if ((handles.snrP(i) >= 10))
                   SNRweightP(i) = 1;
                 elseif ((handles.snrP(i) < 10) && (handles.snrP(i) > 1))
	           % Linear function between two extreme values; NOTE: if change the weighting limits
                   % will need to change this equation too 
                   SNRweightP(i) = (0.1056*handles.snrP(i)) - 0.0556;       
                 elseif ((handles.snrP(i) <= 1))
                   SNRweightP(i) = 0.05;
                 end
               end
             end 

             % Add the accumlated time shift to the full (not cut) P record (using a copy) and make a
             % new cut dataset with the new t1 time
             t1(i) = getheader(handles.dataPcopy(i),'t1');
             newT1(i) = t1(i) + handles.sumShiftP(i);
             handles.dataPcopy(i) = changeheader(handles.dataPcopy(i),'t1',newT1(i));
             handles.dataPcut(i) = cut(handles.dataPcopy(i),'t1',cutMin,'t1',cutMax);
      
             % Before normalize waveform amps, save their max amplitude
             maxCheck = 0;
             for (j = 1:length(handles.dataPcut(i).dep))
               if ((j >= searchWinStart) && (j <= searchWinEnd) && (handles.dataPcut(i).dep(j) > maxCheck))
                 maxCheck = handles.dataPcut(i).dep(j); 
               end
             end
             saveMaxAmpP(i) = maxCheck;

             % Again, normalize amplitudes for cut P wave records so max amplitude = 1
             maxampP = max(handles.dataPcut(i).dep);
             handles.dataPcut(i).dep = (handles.dataPcut(i).dep) / maxampP;
          end
    
          % Move the current weighted stack to be the prior stack (for the next iteration)
          handles = rmfield(handles,'oldPstack');
          handles.oldPstack = handles.weightPstack;
        end
      else
        % Include a check so don't get stuck in an infinite loop; if the weighting/stacking scheme above actually makes the
        % correlation/comparison worse, throw it out of the loop (set compare to a big number to do this).
        compare = 99999;
      end
    end

    % At this point, we have adjusted the P stack as best as we can and we want to use the "old" stack as our final weighted
    % stack (because the most recent new one didn't update anything further
    handles = rmfield(handles,'weightPstack');
    handles.weightPstack = handles.oldPstack;
               
    % Determine average SNR, CCC, and GCARC distance for P-stack
    clearvars snrSum cccSum distSum count snrAvg cccAvg distAvg numSta
    snrSum = 0;
    cccSum = 0;
    distSum = 0;
    count = 0;
    for i=1:length(handles.dataPcut) 
       if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0)
         snrSum = snrSum + handles.snrP(i);
         cccSum = cccSum + handles.coefP(i);
         distSum = distSum + (getheader(handles.dataPcut(i),'GCARC'));
         count = count + 1;
       end 
    end
    snrAvg = strcat('Avg SNR:',32,(num2str(snrSum / count)));
    cccAvg = strcat('Avg CCC:',32,(num2str(cccSum / count)));
    distAvg = strcat('Avg GCARC:',32,(num2str(distSum / count)));
    numSta = strcat('Numb Stations:',32,(num2str(count)));
    
    % Compute P stack standard deviation
    % First, populate a matrix with each column being a weighted P waveform (pointsW); will keep a version of
    % the "non weighted" records too.  This doesn't include any records with SNRweight = 0 though.
    clearvars points pointsW stdDevP lowSnrEvts pupper plower X Y
    for i = (1:length(handles.dataPcut))
      j = 1;
      while (j <= (length(handles.dataPcut(i).dep)))
        if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0)
          points(j,i) = handles.dataPcut(i).dep(j);
          pointsW(j,i) = points(j,i) * (handles.coefP(i)^2) * SNRweightP(i) * handles.stationW(i,1) * handles.evtWP(i,1);
        end
        j = j + 1;
      end
      j = j - 1;
    end
    % Find the number of records that have an SNR value less than 1 (not counting those with SNR weight = 0)
    lowSnrEvts = 0;
    for i = 1:length(handles.snrP)
      if ((handles.snrP(i) < 1) && (SNRweightP(i) ~= 0) && (handles.evtWP(i,1) ~= 0))
        lowSnrEvts = lowSnrEvts + 1;
      end
    end
    % If the number of low SNR records is more than 50% of the total number of stations (disregarding those with
    % SNR = 0), then just compute a normal, non-weighted standard deviation.  Otherwise compute a weighted
    % standard deviation instead.  Doing this because if all records (or many records) have lousy signal and hence
    % low weight, the standard deviation is "articifically" down-weighted too so it looks like the stack has little
    % to no variability.  This is done to get around that issue.   
    if (lowSnrEvts > (count/2))
      for i = 1:length(points)
        stdDevP(i) = std(points(i,:));
      end
    else 
      for i = 1:length(pointsW)
        stdDevP(i) = std(pointsW(i,:));
      end
    end
 
    % Define upper and lower standard deviation limits and associated times
    pupper = handles.weightPstack' + stdDevP;
    plower = handles.weightPstack' - stdDevP;
    X = [timeP, fliplr(timeP)];
    Y = [pupper, fliplr(plower)];
   
    % Plot standard deviation
    cla(handles.axes6,'reset');
    axes(handles.axes6);
    fill(X,Y,[0.65 0.65 0.65],'edgecolor','none');
    hold on

    % Plot weighted P stack on top of standard deviation shading
    plot(timeP,handles.weightPstack,'R','LineWidth',1);
    plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
    xlim([cutMin cutMax]);
    ylim([-1.2 1.2]);
    set(gca,'XMinorTick','on');
    xlabel('Time relative to PREM predicted P time');
    ylabel('Norm. Amplitude');
    rectangle('Position',[-25 0.4 11.5 0.8],'FaceColor','w','EdgeColor','k')
    rectangle('Position',[-25 -1.2 11.5 0.55],'FaceColor','w','EdgeColor','k')
    text(-24.5,1,snrAvg,'color','k','Fontsize',12);
    text(-24.5,0.75,cccAvg,'color','k','Fontsize',12);
    text(-24.5,0.50,distAvg,'color','k','Fontsize',12);
    text(-24.5,-1.1,numSta,'color','k','Fontsize',12);

    % Plot stack with waveforms and report SNR and CCC values for each; if SNR weight = 0, text is shown
    % in gray to indicate that record wasn't used in the stack
    cla(handles.axes1,'reset');
    axes(handles.axes1);
    l = 0;
    m = 1.25;
    clearvars dist snrVal coefVal shiftVal PPtimes evtWeightVal compWeightVal
    for i = 1:length(handles.dataPcut)
       ampP = (handles.dataPcut(i).dep)+l;
       if (SNRweightP(i) == 0 || handles.evtWP(i,1) == 0)
         plot(timeP,ampP,'color',[0.5 0.5 0.5],'LineWidth',1);
       else
         plot(timeP,ampP,'K','LineWidth',1);
         hold on
         plot(timeP,handles.weightPstack+l,'r','LineWidth',1);
       end
       hold on
      
       % Use tauptime to determine PP arrival time along with its reflections/triplications so can add to plots (if they fall in the window)
       PPtimes = tauptime('mod','prem','dep',(getheader(handles.dataP(i),'EVDP')),'sta',[(getheader(handles.dataP(i),'STLA')),(getheader(handles.dataP(i),'STLO'))],'evt',[(getheader(handles.dataP(i),'EVLA')),(getheader(handles.dataP(i),'EVLO'))],'ph','PP,pPP,sPP');
       for r = 1:length(PPtimes)
         PPtt = (((PPtimes(r).time)+(getheader(handles.dataP(i),'O')))-cutMax) - ((getheader(handles.dataP(i),'T1'))-cutMax);
         if ((PPtt >= cutMin) && (PPtt <= cutMax))
           plot([PPtt PPtt],[-1+l 1+l],'color','m','LineWidth',2);
           text(PPtt+0.5,-1.25+l,(PPtimes(r).phase),'color','m','Fontsize',12);
         end
       end
 
       % Add pP and sP markers to records too (provided they fall in the window)
       if ((pickTimes(i,1) >= cutMin) && (pickTimes(i,1) <= cutMax))
         plot([pickTimes(i,1) pickTimes(i,1)],[-1+l 1+l],'color','m','LineWidth',2);
         text(pickTimes(i,1)+0.5,-1.25+l,'pP','color','m','Fontsize',12);
       end
       if ((pickTimes(i,2) >= cutMin) && (pickTimes(i,2) <= cutMax))
         plot([pickTimes(i,2) pickTimes(i,2)],[-1+l 1+l],'color','m','LineWidth',2);
         text(pickTimes(i,2)+0.5,-1.25+l,'sP','color','m','Fontsize',12);
       end 

       dist = strcat('GCARC:',32,(num2str(getheader(handles.dataPcut(i),'GCARC'))));
       snrVal = strcat('SNR:',32,(num2str(handles.snrP(i))));
       coefVal = strcat('CCC:',32,(num2str(handles.coefP(i))));
       shiftVal = strcat('Tot Time Shift:',32,(num2str(handles.sumShiftP(i))));
       evtWeightVal = strcat('Evt Weight:',32,(num2str(handles.evtWP(i,1))));
       compWeightVal = strcat('Comp Weight:',32,(num2str((handles.coefP(i)^2) * SNRweightP(i) * handles.evtWP(i,1) * handles.stationW(i,1))));

       if (SNRweightP(i) == 0 || handles.evtWP(i,1) == 0)
         text(-24,m+2,(getheader(handles.dataPcut(i),'KSTNM')),'color',[0.5 0.5 0.5],'Fontsize',12);
         text(-24,m,dist,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax-8),m+2,snrVal,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax-8),m,coefVal,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax+0.5),m+2.5,shiftVal,'color',[0.5 0.5 0.5],'Fontsize',10);
         text((cutMax+0.5),m+1,evtWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
         text((cutMax+0.5),m-0.5,compWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
       else
         text(-24,m+2,(getheader(handles.dataPcut(i),'KSTNM')),'color','k','Fontsize',12);
         text(-24,m,dist,'color','k','Fontsize',12);
         text((cutMax-8),m+2,snrVal,'color','k','Fontsize',12);
         text((cutMax-8),m,coefVal,'color','k','Fontsize',12);
         text((cutMax+0.5),m+2.5,shiftVal,'color','k','Fontsize',10);
         text((cutMax+0.5),m+1,evtWeightVal,'color','k','Fontsize',10);
         text((cutMax+0.5),m-0.5,compWeightVal,'color','k','Fontsize',10);
       end
       hold on;
       l = l+5;
       m = m+5;
    end
    plot([0 0],[-2 (length(handles.dataPcut)*5)],'LineWidth',0.5,'color',[0.5 0.5 0.5]);
    xlim([cutMin cutMax]);
    set(gca,'XMinorTick','on');
    xlabel('Time relative to PREM predicted P time');
    ylim([-2 (length(handles.dataPcut)*5)]);
    zoom on;	% Allow zooming on seismograms
      
    % Ask user if want to remove any P records and restack (note: alignment doesn't need to be repeated)
    clearvars prompt
    prompt = inputdlg({'Are there any stations you want to remove from the P stack (y or n)?'});
    % If they do, present with list of stations (in same order as plotted) for them to select which ones to remove
    if (strcmp(prompt,'y'))
      i = length(handles.dataPcut);
      j = 1;
      while (i >= 1)
        stationList(j,1) = getheader(handles.dataPcut(i),'KSTNM');
        i = i - 1;
        j = j + 1;
      end
      [selectedSta] = listdlg('PromptString',{'Select Stations to Remove'},'SelectionMode','multiple','ListString',stationList,'ListSize',[160 200]);
      selectedSta = (length(handles.dataPcut)+1) - selectedSta;

      % The "selectedSta" values denote the indices of the records the user wants to remove; set corresponding SNR weight = 0 
      for i = 1:length(handles.dataPcut)
        nix = 0;
        for j = 1:length(selectedSta)
          if (i == selectedSta(j))
            nix = 1;
          end
        end
        if (nix == 1)
          SNRweightP(i) = 0;
        end
      end

      % Clear out old values 
      handles = rmfield(handles,'weightPstack');
      clearvars count weightP;

      % Make a new weighted stack
      handles.weightPstack = 0;
      count = 0;
      for i = 1:length(handles.dataPcut)
        weightP = handles.dataPcut(i).dep * SNRweightP(i) * (handles.coefP(i)^2) * handles.stationW(i,1) * handles.evtWP(i,1);
        if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0)
          handles.weightPstack = handles.weightPstack + weightP;
          count = count + 1;
        end
      end
      handles.weightPstack = handles.weightPstack / count;
      handles.weightPstack = (handles.weightPstack) / (max(handles.weightPstack));

      % Recompute CCC values (between new stack and individual records); use same approach as before
      % First cut/taper weighted P stack
      clearvars searchWinStart searchWinEnd maxampP posA startWin endWin extractPstack taperedPstack ZF ZB TPstack
      searchWinStart = ((length(handles.weightPstack))/2)-limitTimeP;
      searchWinEnd = ((length(handles.weightPstack))/2)+limitTimeP;
      maxampP = 0;
      for (i = 1:length(handles.weightPstack))
        if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPstack(i,1) > maxampP))
          maxampP = handles.weightPstack(i,1);
          posA = i;
        end
      end
      startWin = posA - windowSizeP;
      endWin = posA + windowSizeP;
      extractPstack = handles.weightPstack(startWin:endWin);
      taperedPstack = (extractPstack .* tukeywin((length(extractPstack)),0.25));
      ZF = (zeros([1,(startWin-1)]))';
      ZB = (zeros([1,(length(handles.weightPstack)-endWin)]))';
      TPstack = [ZF;taperedPstack;ZB];

      % Then cut/taper individual P records the same as the stack (use same startWin and endWin); correlate with cut/tapered stack
      % to get new CCC value
      clearvars extractPsignal taperPsignal ZF ZB TPsignal acorr
      for i = 1:length(handles.dataPcut)
        extractPsignal = handles.dataPcut(i).dep(startWin:endWin);
        taperPsignal = (extractPsignal .* tukeywin((length(extractPsignal)),0.25));
        ZF = (zeros([1,(startWin-1)]))';
        ZB = (zeros([1,(length(handles.dataPcut(i).dep)-endWin)]))';
        TPsignal = [ZF;taperPsignal;ZB];
        acorr = xcorr(TPsignal,TPstack,'coeff');
        handles.coefP(i) = max(abs(acorr));
      end

      % Recompute average SNR, CCC, and GCARC distance for P-stack
      clearvars snrSum cccSum distSum count snrAvg cccAvg distAvg numSta
      snrSum = 0;
      cccSum = 0;
      distSum = 0;
      count = 0;
      for i=1:length(handles.dataPcut)
         if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0)
           snrSum = snrSum + handles.snrP(i);
           cccSum = cccSum + handles.coefP(i);
           distSum = distSum + (getheader(handles.dataPcut(i),'GCARC'));
           count = count + 1;
         end 
      end
      snrAvg = strcat('Avg SNR:',32,(num2str(snrSum / count)));
      cccAvg = strcat('Avg CCC:',32,(num2str(cccSum / count)));
      distAvg = strcat('Avg GCARC:',32,(num2str(distSum / count)));
      numSta = strcat('Numb Stations:',32,(num2str(count)));
    
      % Recompute P stack standard deviation
      % First, populate a matrix with each column being a weighted P waveform (pointsW); will keep a version of
      % the "non weighted" records too.  This doesn't include any records with SNRweight = 0 though.
      clearvars points pointsW stdDevP lowSnrEvts pupper plower X Y
      for i = (1:length(handles.dataPcut))
        j = 1;
        while (j <= (length(handles.dataPcut(i).dep)))
          if (SNRweightP(i) ~= 0 && handles.evtWP(i,1) ~= 0)
            points(j,i) = handles.dataPcut(i).dep(j);
            pointsW(j,i) = points(j,i) * (handles.coefP(i)^2) * SNRweightP(i) * handles.stationW(i,1) * handles.evtWP(i,1);
          end
          j = j + 1;
        end
        j = j - 1;
      end
      % Find the number of records that have an SNR value less than 1 (not counting those with SNR weight = 0)
      lowSnrEvts = 0;
      for i = 1:length(handles.snrP)
        if ((handles.snrP(i) < 1) && (SNRweightP(i) ~= 0) && (handles.evtWP(i,1) ~= 0))
          lowSnrEvts = lowSnrEvts + 1;
        end
      end
      % If the number of low SNR records is more than 50% of the total number of stations (disregarding those with
      % SNR = 0), then just compute a normal, non-weighted standard deviation.  Otherwise compute a weighted
      % standard deviation instead.  Doing this because if all records (or many records) have lousy signal and hence
      % low weight, the standard deviation is "articifically" down-weighted too so it looks like the stack has little
      % to no variability.  This is done to get around that issue.   
      if (lowSnrEvts > (count/2))
        for i = 1:length(points)
          stdDevP(i) = std(points(i,:));
        end
      else
        for i = 1:length(pointsW)
          stdDevP(i) = std(pointsW(i,:));
        end
      end
 
      % Define upper and lower standard deviation limits and associated times
      pupper = handles.weightPstack' + stdDevP;
      plower = handles.weightPstack' - stdDevP;
      X = [timeP, fliplr(timeP)];
      Y = [pupper, fliplr(plower)];
   
      % Plot standard deviation
      cla(handles.axes6,'reset');
      axes(handles.axes6);
      fill(X,Y,[0.65 0.65 0.65],'edgecolor','none');
      hold on

      % Plot weighted P stack on top of standard deviation shading
      plot(timeP,handles.weightPstack,'R','LineWidth',1);
      plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
      xlim([cutMin cutMax]);
      ylim([-1.2 1.2]);
      xlabel('Time relative to PREM predicted P time');
      ylabel('Norm. Amplitude');
      rectangle('Position',[-25 0.4 11.5 0.8],'FaceColor','w','EdgeColor','k')
      rectangle('Position',[-25 -1.2 11.5 0.55],'FaceColor','w','EdgeColor','k')
      text(-24.5,1,snrAvg,'color','k','Fontsize',12);
      text(-24.5,0.75,cccAvg,'color','k','Fontsize',12);
      text(-24.5,0.50,distAvg,'color','k','Fontsize',12);
      text(-24.5,-1.1,numSta,'color','k','Fontsize',12);

      % Replot individual records (taking out those that we removed)
      cla(handles.axes1,'reset');
      axes(handles.axes1);
      clearvars ampP
      l = 0;
      m = 1.25;
      clearvars PPtimes 
      for i = 1:length(handles.dataPcut)
        ampP = (handles.dataPcut(i).dep)+l;
        if (SNRweightP(i) == 0 || handles.evtWP(i,1) == 0)
          plot(timeP,ampP,'color',[0.5 0.5 0.5],'LineWidth',1);
        else
          plot(timeP,ampP,'K','LineWidth',1);
          hold on
          plot(timeP,handles.weightPstack+l,'r','LineWidth',1);
        end
        hold on
       
        % Use tauptime to determine PP arrival time along with its reflections/triplications so can add to plots (if they fall in the window)
        PPtimes = tauptime('mod','prem','dep',(getheader(handles.dataP(i),'EVDP')),'sta',[(getheader(handles.dataP(i),'STLA')),(getheader(handles.dataP(i),'STLO'))],'evt',[(getheader(handles.dataP(i),'EVLA')),(getheader(handles.dataP(i),'EVLO'))],'ph','PP,pPP,sPP');
        for r = 1:length(PPtimes)
          PPtt = (((PPtimes(r).time)+(getheader(handles.dataP(i),'O')))-cutMax) - ((getheader(handles.dataP(i),'T1'))-cutMax);
          if ((PPtt >= cutMin) && (PPtt <= cutMax))
            plot([PPtt PPtt],[-1+l 1+l],'color','m','LineWidth',2);
            text(PPtt+0.5,-1.25+l,(PPtimes(r).phase),'color','m','Fontsize',12);
          end
        end
 
        % Add pP and sP markers to records too (provided they fall in the window)
        if ((pickTimes(i,1) >= cutMin) && (pickTimes(i,1) <= cutMax))
          plot([pickTimes(i,1) pickTimes(i,1)],[-1+l 1+l],'color','m','LineWidth',2);
          text(pickTimes(i,1)+0.5,-1.25+l,'pP','color','m','Fontsize',12);
        end
        if ((pickTimes(i,2) >= cutMin) && (pickTimes(i,2) <= cutMax))
          plot([pickTimes(i,2) pickTimes(i,2)],[-1+l 1+l],'color','m','LineWidth',2);
          text(pickTimes(i,2)+0.5,-1.25+l,'sP','color','m','Fontsize',12);
        end 

        dist = strcat('GCARC:',32,(num2str(getheader(handles.dataPcut(i),'GCARC'))));
        snrVal = strcat('SNR:',32,(num2str(handles.snrP(i))));
        coefVal = strcat('CCC:',32,(num2str(handles.coefP(i))));
        shiftVal = strcat('Tot Time Shift:',32,(num2str(handles.sumShiftP(i))));
        evtWeightVal = strcat('Evt Weight:',32,(num2str(handles.evtWP(i,1))));
        compWeightVal = strcat('Comp Weight:',32,(num2str((handles.coefP(i)^2) * SNRweightP(i) * handles.evtWP(i,1) * handles.stationW(i,1))));
       
        if (SNRweightP(i) == 0 || handles.evtWP(i,1) == 0)
          text(-24,m+2,(getheader(handles.dataPcut(i),'KSTNM')),'color',[0.5 0.5 0.5],'Fontsize',12);
          text(-24,m,dist,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax-8),m+2,snrVal,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax-8),m,coefVal,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax+0.5),m+2.5,shiftVal,'color',[0.5 0.5 0.5],'Fontsize',10);
          text((cutMax+0.5),m+1,evtWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
          text((cutMax+0.5),m-0.5,compWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
        else
          text(-24,m+2,(getheader(handles.dataPcut(i),'KSTNM')),'color','k','Fontsize',12);
          text(-24,m,dist,'color','k','Fontsize',12);
          text((cutMax-8),m+2,snrVal,'color','k','Fontsize',12);
          text((cutMax-8),m,coefVal,'color','k','Fontsize',12);
          text((cutMax+0.5),m+2.5,shiftVal,'color','k','Fontsize',10);
          text((cutMax+0.5),m+1,evtWeightVal,'color','k','Fontsize',10);
          text((cutMax+0.5),m-0.5,compWeightVal,'color','k','Fontsize',10);
        end
        l = l+5;
        m = m+5;
      end
       plot([0 0],[-2 (length(handles.dataPcut)*5)],'LineWidth',0.5,'color',[0.5 0.5 0.5]);
       xlim([cutMin cutMax]);
       set(gca,'XMinorTick','on');
       xlabel('Time relative to PREM predicted P time');
       ylim([-2 (length(handles.dataPcut)*5)]);
       zoom on;
    end

    %% Now - move onto PcP data 

    % First, apply the total time shifts computed for the P waveforms to the corresponding PcP waveforms 
    clearvars t4 newT4 maxampPcP
    handles.dataPcP = handles.dataPcPorig;
    for i = (1:length(handles.dataPcP))
      t4(i) = getheader(handles.dataPcP(i),'t4');
      newT4(i) = t4(i) + handles.sumShiftP(i); 
      handles.dataPcP(i) = changeheader(handles.dataPcP(i),'t4',newT4(i));
      handles.dataPcPcut(i) = cut(handles.dataPcP(i),'t4',cutMin,'t4',cutMax);
      maxampPcP = max(handles.dataPcPcut(i).dep);
      handles.dataPcPcut(i).dep = (handles.dataPcPcut(i).dep) / maxampPcP; 
    
      % Also, check to see if the corresponding P wave was "kicked out" of the analysis 
      % (i.e., its time shift was too big).  If so, give its PcP signal an SNR weight = 0
      % so that it is also kicked out
      if (SNRweightP(i) == 0)
        SNRweightPcP(i) = 0;
      end
    end
    
    % Make linear PcP-wave stack (no weighting)
    clearvars pcpStack0 count
    pcpStack0 = 0;
    count = 0;
    for i = 1:length(handles.dataPcPcut)
       if (SNRweightPcP(i) ~= 0)
         pcpStack0 = pcpStack0 + handles.dataPcPcut(i).dep;
         count = count + 1;
       end
    end
    handles.oldPcPstack = pcpStack0 / count;
    handles.oldPcPstack = handles.oldPcPstack / (max(handles.oldPcPstack));

    % Use the same approach as for P above: assume PcP arrives within +/- limitTime of the predicted arrival time
    % (zero on the plots).  Find max amplitude within that window; extract portion of PcP stack +/- windowSizeT
    % around that max amplitude; taper and pad with zeros so length is same as stack
    clearvars searchWinStart searchWinEnd maxampPcP posB startWin endWin extractPcPstack taperedPcPstack ZF ZB TPcPstack
    searchWinStart = ((length(handles.oldPcPstack))/2)-limitTimeP;
    searchWinEnd = ((length(handles.oldPcPstack))/2)+limitTimeP;
    maxampPcP = 0;
    for (i = 1:length(handles.oldPcPstack))
      if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.oldPcPstack(i,1) > maxampPcP))
        maxampPcP = handles.oldPcPstack(i,1);
        posB = i;
      end
    end
    startWin = posB - windowSizeP;
    endWin = posB + windowSizeP;
    extractPcPstack = handles.oldPcPstack(startWin:endWin);
    taperedPcPstack = (extractPcPstack .* tukeywin((length(extractPcPstack)),0.25));
    ZF = (zeros([1,(startWin-1)]))';
    ZB = (zeros([1,(length(handles.oldPcPstack)-endWin)]))';
    TPcPstack = [ZF;taperedPcPstack;ZB];

    % Cut and taper the individual PcP records the same as the stack (use same startWin and endWin)
    clearvars extractPcPsignal taperPcPsignal ZF ZB TPcPsignal acor lag I acorr t4 newT4 maxampPcP 
    handles.dataPcPcopy = handles.dataPcP;
    for i = 1:length(handles.dataPcPcut)
      extractPcPsignal = handles.dataPcPcut(i).dep(startWin:endWin);
      taperPcPsignal = (extractPcPsignal .* tukeywin((length(extractPcPsignal)),0.25));
      ZF = (zeros([1,(startWin-1)]))';
      ZB = (zeros([1,(length(handles.dataPcPcut(i).dep)-endWin)]))';
      TPcPsignal = [ZF;taperPcPsignal;ZB];

      % Cross-correlate the extracted/tapered PcP stack with each extracted/tapered individual PcP record
      % to determine a time shift and the associated CCC values
      [acor,lag] = xcorr(TPcPsignal,TPcPstack);
      [~,I] = max(acor);
      lagDiff(i) = lag(I);
      handles.timeShiftPcP(i) = lagDiff(i) / sampRate;
      acorr = xcorr(TPcPsignal,TPcPstack,'coeff');
      handles.coefPcP(i) = max(abs(acorr));
      handles.sumShiftPcP(i) = handles.sumShiftPcP(i) + handles.timeShiftPcP(i);
    
      % Apply the correlated time shift to adjust the t4 marker, capping the shift at +/- maxShift
      if ((abs(handles.sumShiftPcP(i))) >= maxShift)
        if (handles.sumShiftPcP(i) < 0)
          handles.sumShiftPcP(i) = (-1 * maxShift);
        else
          handles.sumShiftPcP(i) = maxShift;
        end
      end

      t4(i) = getheader(handles.dataPcPcopy(i),'t4');
      newT4(i) = t4(i) + handles.sumShiftPcP(i);
      handles.dataPcPcopy(i) = changeheader(handles.dataPcPcopy(i),'t4',newT4(i));
      handles.dataPcPcut(i) = cut(handles.dataPcPcopy(i),'t4',cutMin,'t4',cutMax);

      % Normalize waveform amps before moving forward
      maxampPcP = max(handles.dataPcPcut(i).dep);
      handles.dataPcPcut(i).dep = (handles.dataPcPcut(i).dep) / maxampPcP;
    end

    % Setup loop to determine SNRs, which will be combined with CCC values to weight each waveform.  The weighted waveforms are 
    % combined into a weighted (normalized) stack.  The loop will continue as long as the new stacks are getting better and until 
    % they are only minimally different from the prior one. 
    clearvars compare compareOld 
    compare = 0.01;
    compareOld = 0;

    while ((compare <= 0.995)) 
      if (compareOld < compare)
        clearvars count 
        handles.weightPcPstack = 0;
        count = 0;
        
        for i = 1:length(handles.dataPcPcut)
           clearvars signaPcPsamp noisePcPsamp signalPcP noisePcP noisePsamp noiseP acor weightPcP
           
           % Compute the signal-to-noise ratio for each individual waveform
           signalPcPsamp = cut(handles.dataPcPcut(i),'t4',signalStart,'t4',signalEnd);
           noisePcPsamp = cut(handles.dataPcPcut(i),'t4',noiseStart,'t4',noiseEnd);
           signalPcP = mean(abs(signalPcPsamp.dep));
           noisePcP = mean(abs(noisePcPsamp.dep));
           handles.snrPcP(i) = signalPcP / noisePcP;
          
           % Regardless of how noise defined above, we are also always going to compute a SNR using the
           % PcP signal and the P noise (to have for comparison).  Signal is same as above but ensure
           % noise is defined correctly here
           noisePsamp = cut(handles.dataPcut(i),'t1',noiseStart,'t1',noiseEnd);
           noiseP = mean(abs(noisePsamp.dep));
           handles.snrPcP2P(i) = signalPcP / noiseP;

           % Setup ramp function for SNR weighting
           if ((handles.snrPcP(i) >= 2) && (SNRweightPcP(i) ~= 0))
             SNRweightPcP(i) = 1;
           elseif ((handles.snrPcP(i) < 2) && (handles.snrPcP(i) > 0.5) && (SNRweightPcP(i) ~= 0))
	     % Linear function between two extreme values; NOTE: if change weighting limits, will 
             % need to change this equation too
             SNRweightPcP(i) = (0.6333*handles.snrPcP(i)) - 0.2667;       
           elseif ((handles.snrPcP(i) <= 0.5) && (SNRweightPcP(i) ~= 0)) 
             SNRweightPcP(i) = 0.05;
           end
     
           weightPcP = handles.dataPcPcut(i).dep * SNRweightPcP(i) * (handles.coefPcP(i)^2) * handles.stationW(i,1) * (handles.evtWP(i,1) * handles.evtWPcP(i,1));
           if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
             handles.weightPcPstack = handles.weightPcPstack + weightPcP;
             count = count + 1;     
           end
        end
        
        % Compute weighted PcP stack; normalize amplitudes
        handles.weightPcPstack = handles.weightPcPstack / count;    
        handles.weightPcPstack = handles.weightPcPstack / (max(handles.weightPcPstack));

        % Compare the newest weighted stack to the previous stack
        acor = xcorr(handles.oldPcPstack,handles.weightPcPstack,'coeff');
        compareOld = compare;
        compare = max(acor);

        % Check similarity between stacks to see if need to keep going
        if (compare <= 0.995)
          % Like above, find max amplitude within specifed time window, extract signal, taper, and use this for the next cross-correlation
          clearvars searchWinStart searchWinEnd maxampPcP posB startWin endWin extractPcPstack taperedPcPstack ZF ZB TPcPstack
          searchWinStart = ((length(handles.weightPcPstack))/2)-limitTimeP;
          searchWinEnd = ((length(handles.weightPcPstack))/2)+limitTimeP;
          maxampPcP = 0;
          for (i = 1:length(handles.weightPcPstack))
            if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPcPstack(i,1) > maxampPcP))
              maxampPcP = handles.weightPcPstack(i,1);
              posB = i;
            end
          end
          startWin = posB - windowSizeP;
          endWin = posB + windowSizeP;
          extractPcPstack = handles.weightPcPstack(startWin:endWin);
          taperedPcPstack = (extractPcPstack .* tukeywin((length(extractPcPstack)),0.25));
          ZF = (zeros([1,(startWin-1)]))';
          ZB = (zeros([1,(length(handles.weightPcPstack)-endWin)]))';
          TPcPstack = [ZF;taperedPcPstack;ZB];

          % Cut and taper the individual PcP records the same as the stack (use same startWin and endWin)
          clearvars extractPcPsignal taperPcPsignal ZF ZB TPcPsignal acor lag I acorr t4 newT4 maxampPcP 
          handles = rmfield(handles,'dataPcPcopy');
          handles.dataPcPcopy = handles.dataPcP;
          for i = 1:length(handles.dataPcPcut)
            extractPcPsignal = handles.dataPcPcut(i).dep(startWin:endWin);
            taperPcPsignal = (extractPcPsignal .* tukeywin((length(extractPcPsignal)),0.25));
            ZF = (zeros([1,(startWin-1)]))';
            ZB = (zeros([1,(length(handles.dataPcPcut(i).dep)-endWin)]))';
            TPcPsignal = [ZF;taperPcPsignal;ZB];

            % Cross-correlate the extracted/tapered PcP stack with each extracted/tapered individual PcP record
            % to determine a time shift and the associated CCC values
            [acor,lag] = xcorr(TPcPsignal,TPcPstack);
            [~,I] = max(acor);
            lagDiff(i) = lag(I);
            handles.timeShiftPcP(i) = lagDiff(i) / sampRate;
            acorr = xcorr(TPcPsignal,TPcPstack,'coeff');
            handles.coefPcP(i) = max(abs(acorr));
            
            % Add this new time shift to the shift sum
            handles.sumShiftPcP(i) = handles.sumShiftPcP(i) + handles.timeShiftPcP(i);
    
            % Apply the correlated time shift to adjust the t4 marker, capping the shift at +/- maxShift
            if ((abs(handles.sumShiftPcP(i))) >= maxShift)
              if (handles.sumShiftPcP(i) < 0)
                handles.sumShiftPcP(i) = (-1 * maxShift);
              else
                handles.sumShiftPcP(i) = maxShift;
              end
            end

            t4(i) = getheader(handles.dataPcPcopy(i),'t4');
            newT4(i) = t4(i) + handles.sumShiftPcP(i);
            handles.dataPcPcopy(i) = changeheader(handles.dataPcPcopy(i),'t4',newT4(i));
            handles.dataPcPcut(i) = cut(handles.dataPcPcopy(i),'t4',cutMin,'t4',cutMax);

            % Normalize waveform amps before moving forward
            maxampPcP = max(handles.dataPcPcut(i).dep);
            handles.dataPcPcut(i).dep = (handles.dataPcPcut(i).dep) / maxampPcP;
          end

          % Move the current weighted stack to be the prior stack (for the next iteration)
          handles = rmfield(handles,'oldPcPstack');
          handles.oldPcPstack = handles.weightPcPstack;
        end
      else
        % Include check that the weighting/shifting process is actually making the correlation better; if not, set compare to
        % a large number and throw out of the loop.  This can happen with bad data and from enforcing the maximum time shift.
        compare = 99999;
      end 
    end

    % At this point, we have adjusted the PcP stack as best as we can and we want to use the "old" stack as our final weighted 
    % stack (because the most recent new one didn't update anything further
    handles = rmfield(handles,'weightPcPstack');
    handles.weightPcPstack = handles.oldPcPstack;
    
    % Determine average SNR, CCC, and GCARC distance for PcP-stack
    clearvars snrSum cccSum distSum count snrAvg cccAvg distAvg numSta
    snrSum = 0;
    cccSum = 0;
    distSum = 0;
    count = 0;
    for i=1:length(handles.dataPcPcut) 
       if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
         snrSum = snrSum + handles.snrPcP(i);
         cccSum = cccSum + handles.coefPcP(i);
         distSum = distSum + (getheader(handles.dataPcPcut(i),'GCARC'));
         count = count + 1;
       end
    end
    snrAvg = strcat('Avg SNR:',32,(num2str(snrSum / count)));
    cccAvg = strcat('Avg CCC:',32,(num2str(cccSum / count)));
    distAvg = strcat('Avg GCARC:',32,(num2str(distSum / count)));
    numSta = strcat('Numb Stations:',32,(num2str(count)));
    
    % Compute PcP stack standard deviation
    % First, populate a matrix with each column being a weighted PcP waveform (pointsW); will keep a version of
    % the "non weighted" records too.  This doesn't include any records with SNRweight = 0 though.
    clearvars points pointsW stdDevPcP lowSnrEvts pupper plower X Y
    for i = (1:length(handles.dataPcPcut))
      j = 1;
      while (j <= (length(handles.dataPcPcut(i).dep)))
        if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
          points(j,i) = handles.dataPcPcut(i).dep(j);
          pointsW(j,i) = points(j,i) * (handles.coefPcP(i)^2) * SNRweightPcP(i) * handles.stationW(i,1) * (handles.evtWP(i,1) * handles.evtWPcP(i,1));
        end 
        j = j + 1;
      end
      j = j - 1;
    end
    % Find the number of records that have an SNR value less than 1 (not counting those with SNR weight = 0)
    lowSnrEvts = 0;
    for i = 1:length(handles.snrPcP)
      if ((handles.snrPcP(i) < 1) && (SNRweightPcP(i) ~= 0) && (handles.evtWP(i,1) ~= 0) && (handles.evtWPcP(i,1) ~= 0))
        lowSnrEvts = lowSnrEvts + 1;
      end
    end
    % If the number of low SNR records is more than 50% of the total number of stations (disregarding those with
    % SNR = 0), then just compute a normal, non-weighted standard deviation.  Otherwise compute a weighted 
    % standard deviation instead.  Doing this because if all records (or many records) have lousy signal and hence
    % low weight, the standard deviation is "artifically" down-weighted too so it looks like the stack has little to
    % no variability.  This is done to get around that issue.
    if (lowSnrEvts > (count/2))
      for i = 1:length(points)
        stdDevPcP(i) = std(points(i,:));
      end
    else
      for i = 1:length(pointsW)
        stdDevPcP(i) = std(pointsW(i,:));
      end
    end

    % Define upper and lower standard deviation limits and associated times
    pupper = handles.weightPcPstack' + stdDevPcP;
    plower = handles.weightPcPstack' - stdDevPcP;
    X = [timePcP, fliplr(timePcP)];
    Y = [pupper, fliplr(plower)];
   
    % Plot standard deviation
    cla(handles.axes7,'reset');
    axes(handles.axes7);
    fill(X,Y,[0.65 0.65 0.65],'edgecolor','none');
    hold on
 
    % Plot weighted PcP stack on top of standard deviation shading
    plot(timePcP,handles.weightPcPstack,'R','LineWidth',1);
    plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
    xlim([cutMin cutMax]);
    ylim([-1.2 1.2]);
    set(gca,'XMinorTick','on');
    xlabel('Time relative to PREM predicted PcP time');
    ylabel('Norm. Amplitude');
    rectangle('Position',[-25 0.4 11.5 0.8],'FaceColor','w','EdgeColor','k')
    rectangle('Position',[-25 -1.2 11.5 0.55],'FaceColor','w','EdgeColor','k')
    text(-24.5,1,snrAvg,'color','k','Fontsize',12);
    text(-24.5,0.75,cccAvg,'color','k','Fontsize',12);
    text(-24.5,0.5,distAvg,'color','k','Fontsize',12);
    text(-24.5,-1.1,numSta,'color','k','Fontsize',12);
 
    % Plot stack with waveforms and report SNR and CCC values; if text is gray, it indicates that record has
    % a SNR weight = 0 and was not included in the stack
    cla(handles.axes2,'reset');
    axes(handles.axes2);
    l = 0;
    m = 1.25;
    clearvars dist snrVal coefVal shiftVal PPtimes evtWeightVal compWeightVal
    for i = 1:length(handles.dataPcPcut)
       ampPcP = (handles.dataPcPcut(i).dep)+l;
       if (SNRweightPcP(i) == 0 || handles.evtWPcP(i,1) == 0 || handles.evtWP(i,1) == 0)
         plot(timePcP,ampPcP,'color',[0.5 0.5 0.5],'LineWidth',1);
       else
         plot(timePcP,ampPcP,'K','LineWidth',1);
         hold on
         plot(timePcP,handles.weightPcPstack+l,'r','LineWidth',1);
       end
       hold on

       % Use tauptime to determine PP arrival time along with its reflection/triplications so can add to plots (if they fall in the window)
       PPtimes = tauptime('mod','prem','dep',(getheader(handles.dataPcPorig(i),'EVDP')),'sta',[(getheader(handles.dataPcPorig(i),'STLA')),(getheader(handles.dataPcPorig(i),'STLO'))],'evt',[(getheader(handles.dataPcPorig(i),'EVLA')),(getheader(handles.dataPcPorig(i),'EVLO'))],'ph','PP,pPP,sPP');
       for r = 1:length(PPtimes)
         PPtt = (((PPtimes(r).time)+(getheader(handles.dataPcPorig(i),'O')))-cutMax) - ((getheader(handles.dataPcPorig(i),'T4'))-cutMax);
         if ((PPtt >= cutMin) && (PPtt <= cutMax))
           plot([PPtt PPtt],[-1+l 1+l],'color','m','LineWidth',2);
           text(PPtt+0.5,-1.25+l,(PPtimes(r).phase),'color','m','Fontsize',12);
         end
       end

       % Add pP and sP markers to records too (provided they fall in the window)
       if ((pickTimes(i,3) >= cutMin) && (pickTimes(i,3) <= cutMax))
         plot([pickTimes(i,3) pickTimes(i,3)],[-1+l 1+l],'color','m','LineWidth',2);
         text(pickTimes(i,3)+0.5,-1.25+l,'pP','color','m','Fontsize',12);
       end
       if ((pickTimes(i,4) >= cutMin) && (pickTimes(i,4) <= cutMax))
         plot([pickTimes(i,4) pickTimes(i,4)],[-1+l 1+l],'color','m','LineWidth',2);
         text(pickTimes(i,4)+0.5,-1.25+l,'sP','color','m','Fontsize',12);
       end
        
       snrVal = strcat('SNR:',32,(num2str(handles.snrPcP(i))));
       dist = strcat('GCARC:',32,(num2str(getheader(handles.dataPcPcut(i),'GCARC'))));
       coefVal = strcat('CCC:',32,(num2str(handles.coefPcP(i))));
       shiftVal = strcat('Tot Time Shift:',32,(num2str(handles.sumShiftPcP(i))));
       evtWeightVal = strcat('Evt Weight:',32,(num2str((handles.evtWP(i,1)*handles.evtWPcP(i,1)))));
       compWeightVal = strcat('Comp Weight:',32,(num2str((handles.coefPcP(i)^2) * SNRweightPcP(i) * (handles.evtWP(i,1) * handles.evtWPcP(i,1)) * handles.stationW(i,1))));

       if (SNRweightPcP(i) == 0 || handles.evtWPcP(i,1) == 0 || handles.evtWP(i,1) == 0)
         text(-24,m+2,(getheader(handles.dataPcPcut(i),'KSTNM')),'color',[0.5 0.5 0.5],'Fontsize',12);
         text(-24,m,dist,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax-9.25),m+2,snrVal,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax-9.25),m-0.25,coefVal,'color',[0.5 0.5 0.5],'Fontsize',12);
         text((cutMax+0.5),m+2.5,shiftVal,'color',[0.5 0.5 0.5],'Fontsize',10);
         text((cutMax+0.5),m+1,evtWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
         text((cutMax+0.5),m-0.5,compWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
       else
         text(-24,m+2,(getheader(handles.dataPcPcut(i),'KSTNM')),'color','k','Fontsize',12);
         text(-24,m,dist,'color','k','Fontsize',12);
         text((cutMax-9.25),m+2,snrVal,'color','k','Fontsize',12);
         text((cutMax-9.25),m-0.25,coefVal,'color','k','Fontsize',12);
         text((cutMax+0.5),m+2.5,shiftVal,'color','k','Fontsize',10);
         text((cutMax+0.5),m+1,evtWeightVal,'color','k','Fontsize',10);
         text((cutMax+0.5),m-0.5,compWeightVal,'color','k','Fontsize',10);
       end
       l = l+5;
       m = m+5;
    end
    plot([0 0],[-2 (length(handles.dataPcPcut)*5)],'LineWidth',0.5,'color',[0.5 0.5 0.5]);
    xlim([cutMin cutMax]);
    set(gca,'XMinorTick','on');
    xlabel('Time relative to PREM predicted PcP time');
    ylim([-2 (length(handles.dataPcPcut)*5)]);
    zoom on;	% Allow zooming on seismograms

    % Like with the P records, ask user if want to remove any PcP records and restack (note: alignment doesn't need to be repeated)
    clearvars prompt
    prompt = inputdlg({'Are there any stations you want to remove from the PcP stack (y or no)?'});
    % If they do, present with list of stations (in same order as plotted) for them to select which ones to remove
    if (strcmp(prompt,'y'))
      clearvars stationList selectedSta
      i = length(handles.dataPcPcut);
      j = 1;
      while (i >= 1)
        stationList(j,1) = getheader(handles.dataPcPcut(i),'KSTNM');
        i = i - 1;
        j = j + 1;
      end
      [selectedSta] = listdlg('PromptString',{'Select Stations to Remove'},'SelectionMode','multiple','ListString',stationList,'ListSize',[160 200]);
      selectedSta = (length(handles.dataPcPcut)+1) - selectedSta;

      % The "selectedSta" values denote the indices of the records the user wants to remove; set their corresponding SNR weight = 0 
      for i = 1:length(handles.dataPcPcut)
        nix = 0;
        for j = 1:length(selectedSta)
          if (i == selectedSta(j))
            nix = 1;
          end
        end
        if (nix == 1)
          SNRweightPcP(i) = 0;
        end
      end

      % Clear out old values 
      handles = rmfield(handles,'weightPcPstack');
      clearvars count weightPcP;

      % Make a new weighted stack with the remaining, aligned records
      handles.weightPcPstack = 0;
      count = 0;
      for i = 1:length(handles.dataPcPcut)
        weightPcP = handles.dataPcPcut(i).dep * SNRweightPcP(i) * (handles.coefPcP(i)^2) * handles.stationW(i,1) * handles.evtWP(i,1) * handles.evtWPcP(i,1);
        if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
          handles.weightPcPstack = handles.weightPcPstack + weightPcP;
          count = count + 1;
        end
      end
      handles.weightPcPstack = handles.weightPcPstack / count;
      handles.weightPcPstack = (handles.weightPcPstack) / (max(handles.weightPcPstack));

      % Recompute CCC values (between new stack and individual records); use same approach as before
      % First cut/taper weighted PcP stack
      clearvars searchWinStart searchWinEnd maxampPcP posB startWin endWin extractPcPstack taperedPcPstack ZF ZB TPcPstack
      searchWinStart = ((length(handles.weightPcPstack))/2)-limitTimeP;
      searchWinEnd = ((length(handles.weightPcPstack))/2)+limitTimeP;
      maxampPcP = 0;
      for (i = 1:length(handles.weightPcPstack))
        if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPcPstack(i,1) > maxampPcP))
          maxampPcP = handles.weightPcPstack(i,1);
          posB = i;
        end
      end
      startWin = posB - windowSizeP;
      endWin = posB + windowSizeP;
      extractPcPstack = handles.weightPcPstack(startWin:endWin);
      taperedPcPstack = (extractPcPstack .* tukeywin((length(extractPcPstack)),0.25));
      ZF = (zeros([1,(startWin-1)]))';
      ZB = (zeros([1,(length(handles.weightPcPstack)-endWin)]))';
      TPcPstack = [ZF;taperedPcPstack;ZB];

      % Then cut/taper individual PcP records the same as the stack (use same startWin and endWin); correlated with tapered stack
      clearvars extractPcPsignal taperedPcPsignal ZF ZB TPcPsignal acorr
      for i = 1:length(handles.dataPcPcut)
        extractPcPsignal = handles.dataPcPcut(i).dep(startWin:endWin);
        taperedPcPsignal = (extractPcPsignal .* tukeywin((length(extractPcPsignal)),0.25));
        ZF = (zeros([1,(startWin-1)]))';
        ZB = (zeros([1,(length(handles.dataPcPcut(i).dep)-endWin)]))';
        TPcPsignal = [ZF;taperedPcPsignal;ZB];
        acorr = xcorr(TPcPsignal,TPcPstack,'coeff');
        handles.coefPcP(i) = max(abs(acorr));
      end

      % Recompute average SNR, CCC, and GCARC distance for PcP-stack
      clearvars snrSum cccSum distSum count snrAvg cccAvg distAvg numSta
      snrSum = 0;
      cccSum = 0;
      distSum = 0;
      count = 0;
      for i = 1:length(handles.dataPcPcut) 
         if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
           snrSum = snrSum + handles.snrPcP(i);
           cccSum = cccSum + handles.coefPcP(i);
           distSum = distSum + (getheader(handles.dataPcPcut(i),'GCARC'));
           count = count + 1;
         end 
      end
      snrAvg = strcat('Avg SNR:',32,(num2str(snrSum / count)));
      cccAvg = strcat('Avg CCC:',32,(num2str(cccSum / count)));
      distAvg = strcat('Avg GCARC:',32,(num2str(distSum / count)));
      numSta = strcat('Numb Stations:',32,(num2str(count)));
    
      % Recompute PcP stack standard deviation
      % First, populate a matrix with each column being a weighted PcP waveform (pointsW); will keep a version of
      % the "non weighted" records too.  This doesn't include any records with SNRweight = 0 though.
      clearvars points pointsW stdDevPcP lowSnrEvts pupper plower X Y
      for i = (1:length(handles.dataPcPcut))
        j = 1;
        while (j <= (length(handles.dataPcPcut(i).dep)))
          if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
            points(j,i) = handles.dataPcPcut(i).dep(j);
            pointsW(j,i) = points(j,i) * (handles.coefPcP(i)^2) * SNRweightPcP(i) * handles.stationW(i,1) * (handles.evtWP(i,1) * handles.evtWPcP(i,1));
          end
          j = j + 1;
        end
        j = j - 1;
      end
      % Find the number of records that have an SNR value less than 1 (not counting those with SNR weight = 0)
      lowSnrEvts = 0;
      for i = 1:length(handles.snrPcP)
        if ((handles.snrPcP(i) < 1) && (SNRweightPcP(i) ~= 0) && (handles.evtWP(i,1) ~= 0) && (handles.evtWPcP(i,1) ~= 0))
          lowSnrEvts = lowSnrEvts + 1;
        end
      end
      % If the number of low SNR records is more than 50% of the total number of stations (disregarding those with
      % SNR = 0), then just compute a normal, non-weighted standard deviation.  Otherwise compute a weighted
      % standard deviation instead.  Doing this because if all records (or many records) have lousy signal and hence
      % low weight, the standard deviation is "articifically" down-weighted too so it looks like the stack has little
      % to no variability.  This is done to get around that issue.   
      if (lowSnrEvts > (count/2))
        for i = 1:length(points)
          stdDevPcP(i) = std(points(i,:));
        end
      else 
        for i = 1:length(pointsW)
          stdDevPcP(i) = std(pointsW(i,:));
        end
      end
 
      % Define upper and lower standard deviation limits and associated times
      pupper = handles.weightPcPstack' + stdDevPcP;
      plower = handles.weightPcPstack' - stdDevPcP;
      X = [timePcP, fliplr(timePcP)];
      Y = [pupper, fliplr(plower)];
   
      % Plot standard deviation
      cla(handles.axes7,'reset');
      axes(handles.axes7);
      fill(X,Y,[0.65 0.65 0.65],'edgecolor','none');
      hold on

      % Plot weighted PcP stack on top of standard deviation shading
      plot(timePcP,handles.weightPcPstack,'R','LineWidth',1);
      plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
      xlim([cutMin cutMax]);
      set(gca,'XMinorTick','on');
      ylim([-1.2 1.2]);
      xlabel('Time relative to PREM predicted PcP time');
      ylabel('Norm. Amplitude');
      rectangle('Position',[-25 0.4 11.5 0.8],'FaceColor','w','EdgeColor','k')
      rectangle('Position',[-25 -1.2 11.5 0.55],'FaceColor','w','EdgeColor','k')
      text(-24.5,1,snrAvg,'color','k','Fontsize',12);
      text(-24.5,0.75,cccAvg,'color','k','Fontsize',12);
      text(-24.5,0.50,distAvg,'color','k','Fontsize',12);
      text(-24.5,-1.1,numSta,'color','k','Fontsize',12);

      % Replot individual records (taking out those that we removed)
      cla(handles.axes2,'reset');
      axes(handles.axes2);
      clearvars ampPcP
      l = 0;
      m = 1.25;
      clearvars PPtimes 
      for i = 1:length(handles.dataPcPcut)
        ampPcP = (handles.dataPcPcut(i).dep)+l;
        if (SNRweightPcP(i) == 0 || handles.evtWPcP(i,1) == 0 || handles.evtWP(i,1) == 0)
          plot(timePcP,ampPcP,'color',[0.5 0.5 0.5],'LineWidth',1);
        else
          plot(timePcP,ampPcP,'K','LineWidth',1);
          hold on
          plot(timePcP,handles.weightPcPstack+l,'r','LineWidth',1);
        end
        hold on
       
        % Use tauptime to determine PP arrival time along with its reflection/triplications so can add to plots (if they fall in the window)
        PPtimes = tauptime('mod','prem','dep',(getheader(handles.dataPcPorig(i),'EVDP')),'sta',[(getheader(handles.dataPcPorig(i),'STLA')),(getheader(handles.dataPcPorig(i),'STLO'))],'evt',[(getheader(handles.dataPcPorig(i),'EVLA')),(getheader(handles.dataPcPorig(i),'EVLO'))],'ph','PP,pPP,sPP');
        for r = 1:length(PPtimes)
          PPtt = (((PPtimes(r).time)+(getheader(handles.dataPcPorig(i),'O')))-cutMax) - ((getheader(handles.dataPcPorig(i),'T4'))-cutMax);
          if ((PPtt >= cutMin) && (PPtt <= cutMax))
            plot([PPtt PPtt],[-1+l 1+l],'color','m','LineWidth',2);
            text(PPtt+0.5,-1.25+l,(PPtimes(r).phase),'color','m','Fontsize',12);
          end
        end

        % Add pP and sP markers to records too (provided they fall in the window)
        if ((pickTimes(i,3) >= cutMin) && (pickTimes(i,3) <= cutMax))
          plot([pickTimes(i,3) pickTimes(i,3)],[-1+l 1+l],'color','m','LineWidth',2);
          text(pickTimes(i,3)+0.5,-1.25+l,'pP','color','m','Fontsize',12);
        end
        if ((pickTimes(i,4) >= cutMin) && (pickTimes(i,4) <= cutMax))
          plot([pickTimes(i,4) pickTimes(i,4)],[-1+l 1+l],'color','m','LineWidth',2);
          text(pickTimes(i,4)+0.5,-1.25+l,'sP','color','m','Fontsize',12);
        end
        
        snrVal = strcat('SNR:',32,(num2str(handles.snrPcP(i))));
        coefVal = strcat('CCC:',32,(num2str(handles.coefPcP(i))));
        shiftVal = strcat('Tot Time Shift:',32,(num2str(handles.sumShiftPcP(i))));
        evtWeightVal = strcat('Evt Weight:',32,(num2str((handles.evtWP(i,1)*handles.evtWPcP(i,1)))));
        compWeightVal = strcat('Comp Weight:',32,(num2str((handles.coefPcP(i)^2) * SNRweightPcP(i) * handles.evtWP(i,1) * handles.evtWPcP(i,1) * handles.stationW(i,1))));
       
        if (SNRweightPcP(i) == 0 || handles.evtWPcP(i,1) == 0 || handles.evtWP(i,1) == 0)
          text(-24,m+2,(getheader(handles.dataPcPcut(i),'KSTNM')),'color',[0.5 0.5 0.5],'Fontsize',12);
          text(-24,m,dist,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax-8),m+2,snrVal,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax-8),m,coefVal,'color',[0.5 0.5 0.5],'Fontsize',12);
          text((cutMax+0.5),m+2.5,shiftVal,'color',[0.5 0.5 0.5],'Fontsize',10);
          text((cutMax+0.5),m+1,evtWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
          text((cutMax+0.5),m-0.5,compWeightVal,'color',[0.5 0.5 0.5],'Fontsize',10);
        else
          text(-24,m+2,(getheader(handles.dataPcPcut(i),'KSTNM')),'color','k','Fontsize',12);
          text(-24,m,dist,'color','k','Fontsize',12);
          text((cutMax-8),m+2,snrVal,'color','k','Fontsize',12);
          text((cutMax-8),m,coefVal,'color','k','Fontsize',12);
          text((cutMax+0.5),m+2.5,shiftVal,'color','k','Fontsize',10);
          text((cutMax+0.5),m+1,evtWeightVal,'color','k','Fontsize',10);
          text((cutMax+0.5),m-0.5,compWeightVal,'color','k','Fontsize',10);
        end
        l = l+5;
        m = m+5;
      end
       plot([0 0],[-2 (length(handles.dataPcPcut)*5)],'LineWidth',0.5,'color',[0.5 0.5 0.5]);
       xlim([cutMin cutMax]);
       set(gca,'XMinorTick','on');
       xlabel('Time relative to PREM predicted PcP time');
       ylim([-2 (length(handles.dataPcPcut)*5)]);
       zoom on;
    end

    %% Now, do some P-to-PcP comparison

    % Similar to above, will limit assessment of max P amplitude to a window +/- (limitTime) around the
    % predicted P time (zero on plots).  Find the max amplitude within that window and mark on plot (blue star) 
    clearvars searchWinStart searchWinEnd maxampP sampMaxP timeMaxP 
    searchWinStart = ((length(handles.weightPstack))/2)-limitTimeP;
    searchWinEnd = ((length(handles.weightPstack))/2)+limitTimeP;
    maxampP = 0;
    for (i = 1:length(handles.weightPstack))
      if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPstack(i,1) > maxampP))
 	maxampP = handles.weightPstack(i,1);
        sampMaxP = i; 
        timeMaxP = (sampMaxP / sampRate) - cutMax;
      end
    end
    axes(handles.axes6);
    hold on
    bsP = plot(timeMaxP,maxampP,'b*','MarkerSize',9,'LineWidth',1);
  
    % Similar to above, will limit assessment of max PcP amplitude to a window +/- (limitTime) around the
    % predicted PcP time (zero on plots).  Find the max amplitude within that window and mark on plot (blue star) 
    clearvars searchWinStart searchWinEnd maxampPcP sampMaxPcP timeMaxPcP 
    searchWinStart = ((length(handles.weightPcPstack))/2)-limitTimeP;
    searchWinEnd = ((length(handles.weightPcPstack))/2)+limitTimeP;
    maxampPcP = 0;
    for (i = 1:length(handles.weightPcPstack))
      if ((i >= searchWinStart) && (i <= searchWinEnd) && (handles.weightPcPstack(i,1) > maxampPcP))
 	maxampPcP = handles.weightPcPstack(i,1);
        sampMaxPcP = i; 
        timeMaxPcP = (sampMaxPcP / sampRate) - cutMax;
      end
    end
    axes(handles.axes7);
    hold on
    bsPcP = plot(timeMaxPcP,maxampPcP,'b*','MarkerSize',9,'LineWidth',1);

    % Find time difference between P and PcP arrivals (blue stars marked on stack plots) and time shift PcP stack to align with
    % P; need to pad PcP stack on one end and trim on the other end to maintain correct number of points 
    clearvars timeDiff padding
    timeDiff = timeMaxP - timeMaxPcP;
    if (timeDiff > 0)
       padding = zeros((((abs(timeDiff))*sampRate)),1); 
       handles.weightPcPstackPad = [padding;handles.weightPcPstack];
       handles.weightPcPstackPad = handles.weightPcPstackPad(1:(length(handles.weightPcPstackPad) - (abs(timeDiff)*sampRate)),1);
    elseif (timeDiff < 0)
       padding = zeros((((abs(timeDiff))*sampRate))-1,1); 
       handles.weightPcPstackPad = [handles.weightPcPstack;padding];
       handles.weightPcPstackPad = handles.weightPcPstackPad((abs(timeDiff)*sampRate):length(handles.weightPcPstackPad),1);
    else
       handles.weightPcPstackPad = handles.weightPcPstack;
    end
    cla(handles.axes8,'reset')
    axes(handles.axes8);

    % Plot together so max amplitudes point overlay
    plot(timeP,handles.weightPstack,'k','LineWidth',1);
    hold on
    plot(timeP,handles.weightPcPstackPad,'r','LineWidth',1);
    plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
    xlim([cutMin cutMax]);
    set(gca,'XMinorTick','on');
    ylim([-1.2 1.2]); 
    xlabel('Time relative to PREM predicted P time');
    ylabel('Norm. Amplitude');

    % Ask user if want to change blue star position on P-stack
    clearvars prompt
    prompt = inputdlg({'Do you want to change the marker on the P stack (y or n)?'});
    % If they do, ask them to enter min/max values that bracket position of interest 
    if (strcmp(prompt,'y'))
      while (strcmp(prompt,'y'))
        axes(handles.axes6);
        hold on
        delete(bsP); 
        clearvars B1 B2 maxampP sampMaxP timeMaxP
        prompt = inputdlg({'Enter minimum bracket time:';'Enter maximum bracket time:'});
        B1 = ((str2num(prompt{1,1}))*sampRate)+((length(handles.weightPstack))/2);
        B2 = ((str2num(prompt{2,1}))*sampRate)+((length(handles.weightPstack))/2);
      
        % Find max amplitude within bracketed times provided by user; this is new blue star point
        maxampP = 0;
        for (i = 1:length(handles.weightPstack))
          if ((i >= B1) && (i <= B2) && (handles.weightPstack(i,1) > maxampP))
            maxampP = handles.weightPstack(i,1);
            sampMaxP = i;
            timeMaxP = (sampMaxP / sampRate) - cutMax;
          end
        end
        bsP = plot(timeMaxP,maxampP,'b*','MarkerSize',9,'LineWidth',1);  
        
        % Replot comparison between P and PcP 
        clearvars timeDiff padding
        timeDiff = timeMaxP - timeMaxPcP;
        if (timeDiff > 0)
          padding = zeros((((abs(timeDiff))*sampRate)),1); 
          handles.weightPcPstackPad = [padding;handles.weightPcPstack];
          handles.weightPcPstackPad = handles.weightPcPstackPad(1:(length(handles.weightPcPstackPad) - (abs(timeDiff)*sampRate)),1);
       elseif (timeDiff < 0)
          padding = zeros((((abs(timeDiff))*sampRate))-1,1); 
          handles.weightPcPstackPad = [handles.weightPcPstack;padding];
          handles.weightPcPstackPad = handles.weightPcPstackPad((abs(timeDiff)*sampRate):length(handles.weightPcPstackPad),1);
       else
          handles.weightPcPstackPad = handles.weightPcPstack;
       end
       cla(handles.axes8,'reset')
       axes(handles.axes8);
       plot(timeP,handles.weightPstack,'k','LineWidth',1);
       hold on
       plot(timeP,handles.weightPcPstackPad,'r','LineWidth',1);
       plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
       xlim([cutMin cutMax]);
       set(gca,'XMinorTick','on');
       ylim([-1.2 1.2]); 
       xlabel('Time relative to PREM predicted P time');
       ylabel('Norm. Amplitude');
   
       prompt = inputdlg({'Adjust P marker again? (y or n)?'});
      end
    end
    
    % Do same thing for PcP-stack
    clearvars prompt
    prompt = inputdlg({'Do you want to change the marker on the PcP stack (y or n)?'});
    % If they do, ask them to enter min/max values that bracket position of interest 
    if (strcmp(prompt,'y'))
      while (strcmp(prompt,'y'))
        axes(handles.axes7);
        hold on
        delete(bsPcP);
        clearvars B1 B2 maxampPcP sampMaxPcP timeMaxPcP
        prompt = inputdlg({'Enter minimum bracket time:';'Enter maximum bracket time:'});
        B1 = ((str2num(prompt{1,1}))*sampRate)+((length(handles.weightPcPstack))/2);
        B2 = ((str2num(prompt{2,1}))*sampRate)+((length(handles.weightPcPstack))/2);
      
        % Find max amplitude within bracketed times provided by user; this is new blue star point
        maxampPcP = 0;
        for (i = 1:length(handles.weightPcPstack))
          if ((i >= B1) && (i <= B2) && (handles.weightPcPstack(i,1) > maxampPcP))
            maxampPcP = handles.weightPcPstack(i,1);
            sampMaxPcP = i;
            timeMaxPcP = (sampMaxPcP / sampRate) - cutMax;
          end
        end
        bsPcP = plot(timeMaxPcP,maxampPcP,'b*','MarkerSize',9,'LineWidth',1);  
        
        % Replot comparison between P and PcP 
        clearvars timeDiff padding
        timeDiff = timeMaxP - timeMaxPcP;
        if (timeDiff > 0)
          padding = zeros((((abs(timeDiff))*sampRate)),1); 
          handles.weightPcPstackPad = [padding;handles.weightPcPstack];
          handles.weightPcPstackPad = handles.weightPcPstackPad(1:(length(handles.weightPcPstackPad) - (abs(timeDiff)*sampRate)),1);
       elseif (timeDiff < 0)
          padding = zeros((((abs(timeDiff))*sampRate))-1,1); 
          handles.weightPcPstackPad = [handles.weightPcPstack;padding];
          handles.weightPcPstackPad = handles.weightPcPstackPad((abs(timeDiff)*sampRate):length(handles.weightPcPstackPad),1);
       else
          handles.weightPcPstackPad = handles.weightPcPstack;
       end
       cla(handles.axes8,'reset')
       axes(handles.axes8);
       plot(timeP,handles.weightPstack,'k','LineWidth',1);
       hold on
       plot(timeP,handles.weightPcPstackPad,'r','LineWidth',1);
       plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
       xlim([cutMin cutMax]);
       set(gca,'XMinorTick','on');
       ylim([-1.2 1.2]); 
       xlabel('Time relative to PREM predicted P time');
       ylabel('Norm. Amplitude');
        
       prompt = inputdlg({'Adjust PcP marker again? (y or n)?'});
      end
    end

    % Use blue star to define signal and noise windows for the P stack; compute SNR and add to plot
    clearvars extractPstackSignal extractPstackNoise PstackSNR stackSnrVal STS STE NTS NTE
    extractPstackSignal = handles.weightPstack((sampMaxP+(signalStartStack*sampRate)):(sampMaxP+(signalEndStack*sampRate)));
    if ((sampMaxP+(noiseStartStack*sampRate)) < 0)
      % If noise window extends past the front end of the stack, then just use the first point on the stack as
      % the start of the noise window
      extractPstackNoise =  handles.weightPstack(1:(sampMaxP+(noiseEndStack*sampRate)));
    else
      extractPstackNoise =  handles.weightPstack((sampMaxP+(noiseStartStack*sampRate)):(sampMaxP+(noiseEndStack*sampRate)));
    end
    PstackSNR = (mean(abs(extractPstackSignal))) / (mean(abs(extractPstackNoise)));
    stackSnrVal = strcat('Stack SNR:',32,(num2str(PstackSNR))); 
    axes(handles.axes6);
    text(-24.5,-0.85,stackSnrVal,'color','k','Fontsize',12);

    % Mark signal and noise window limits on stack plot 
    STS = ((sampMaxP+(signalStartStack*sampRate)) / sampRate) - cutMax;
    STE = ((sampMaxP+(signalEndStack*sampRate)) / sampRate) - cutMax;
    NTS = ((sampMaxP+(noiseStartStack*sampRate)) / sampRate) - cutMax;
    NTE = ((sampMaxP+(noiseEndStack*sampRate)) / sampRate) - cutMax;
    plot([STS STS],[-0.25 0.25],'b','LineWidth',1); 
    plot([STE STE],[-0.25 0.25],'b','LineWidth',1);
    plot([NTS NTS],[-0.25 0.25],'b','LineWidth',1);
    plot([NTE NTE],[-0.25 0.25],'b','LineWidth',1); 

    % Use blue star to define signal and noise windows for the PcP stack; compute SNR and add to plot
    clearvars extractPcPstackSignal extractPcPstackNoise PcPstackSNR stackSnrVal STS STE NTS NTE
    extractPcPstackSignal = handles.weightPcPstack((sampMaxPcP+(signalStartStack*sampRate)):(sampMaxPcP+(signalEndStack*sampRate)));
    if ((sampMaxPcP+(noiseStartStack*sampRate)) < 0)
      % If noise window extends past the front end of the stack, then just use the first point on the stack as  
      % the start of the noise window
      extractPcPstackNoise =  handles.weightPcPstack(1:(sampMaxPcP+(noiseEndStack*sampRate)));
    else
      extractPcPstackNoise =  handles.weightPcPstack((sampMaxPcP+(noiseStartStack*sampRate)):(sampMaxPcP+(noiseEndStack*sampRate)));
    end 
    PcPstackSNR = (mean(abs(extractPcPstackSignal))) / (mean(abs(extractPcPstackNoise)));
    stackSnrVal = strcat('Stack SNR:',32,(num2str(PcPstackSNR))); 
    axes(handles.axes7);
    text(-24.5,-0.85,stackSnrVal,'color','k','Fontsize',12);

    % Mark signal and noise window limits on stack plot 
    STS = ((sampMaxPcP+(signalStartStack*sampRate)) / sampRate) - cutMax;
    STE = ((sampMaxPcP+(signalEndStack*sampRate)) / sampRate) - cutMax;
    NTS = ((sampMaxPcP+(noiseStartStack*sampRate)) / sampRate) - cutMax;
    NTE = ((sampMaxPcP+(noiseEndStack*sampRate)) / sampRate) - cutMax;
    plot([STS STS],[-0.25 0.25],'b','LineWidth',1); 
    plot([STE STE],[-0.25 0.25],'b','LineWidth',1);
    plot([NTS NTS],[-0.25 0.25],'b','LineWidth',1);
    plot([NTE NTE],[-0.25 0.25],'b','LineWidth',1); 

    % Find difference between the time-adjusted PcP stack and the P stack ("remainder trace"); note: amps not normalized
    clearvars diffAmp
    diffAmp = handles.weightPcPstackPad - handles.weightPstack;
    cla(handles.axes12,'reset')
    axes(handles.axes12);
    plot(timeP,diffAmp,'b','LineWidth',1);
    hold on
    plot([0 0],[-1.2 1.2],'LineWidth',0.5,'color',[0.5 0.5 0.5])
    xlim([cutMin cutMax]);
    set(gca,'XMinorTick','on');
    ylim([-1.2 1.2]); 
    xlabel('Time relative to PREM predicted P time');
    ylabel('Residual Amplitude');

    % Look at PcP/P amplitude ratios; make copy of original (uncut) PcP waveforms; apply total correlated time shift to
    % shift records; make new cut PcP records that have their respective SNR and CCC weights applied; divide by the max
    % amplitude of the corresponding P record; use these waveforms to make a new weighted PcP stack
    handles = rmfield(handles,'dataPcPcopy');
    clearvars t4 newT4 ARsum count ARstack 
    handles.dataPcPcopy = handles.dataPcPorig;
    ARsum = 0;
    count = 0;

    for i = 1:length(handles.dataPcPcopy)
      t4(i) = getheader(handles.dataPcPcopy(i),'t4');
      newT4(i) = t4(i) + handles.sumShiftP(i) + handles.sumShiftPcP(i);
      handles.dataPcPcopy(i) = changeheader(handles.dataPcPcopy(i),'t4',newT4(i));
      handles.dataPcPcut(i) = cut(handles.dataPcPcopy(i),'t4',cutMin,'t4',cutMax);
      handles.dataPcPcut(i).dep = handles.dataPcPcut(i).dep * (handles.coefPcP(i)^2) * SNRweightPcP(i) * handles.stationW(i,1) * handles.evtWP(i,1) * handles.evtWPcP(i,1); 
      handles.dataPcPcut(i).dep = handles.dataPcPcut(i).dep / saveMaxAmpP(i);
      
      if (SNRweightPcP(i) ~= 0 && handles.evtWP(i,1) ~= 0 && handles.evtWPcP(i,1) ~= 0)
        % Make a weighted stack using all records whose SNR weight is not equal to one
        ARsum = ARsum + handles.dataPcPcut(i).dep;
        count = count + 1;
      end
    end

    % Compute weighted stack with max P amplitude factored in
    ARstack = ARsum / count;

    % Find the max amplitude of the AR stack within the specified search window
    clearvars searchWinStart searchWinEnd maxamp ratioText
    searchWinStart = ((length(ARstack))/2)-limitTimeP;
    searchWinEnd = ((length(ARstack))/2)+limitTimeP;
    maxamp = 0;
    for (i = 1:length(ARstack))
      if ((i >= searchWinStart) && (i <= searchWinEnd) && (ARstack(i,1) > maxamp))
 	maxamp = ARstack(i,1);
      end
    end
    ratioText = strcat('PcP/P amplitude ratio:',32,(num2str(maxamp)));
    axes(handles.axes7);
    text(10,1.35,ratioText,'color','k','FontSize',12);

    % Compute and report the average PcP-P time now that records have been time shifted
    handles = rmfield(handles,'dataPcPcopy');
    handles = rmfield(handles,'dataPcopy');
    clearvars t1 t4 newT1 newT4 sumPhaseTimeDiff
    sumPhaseTimeDiff = 0;
    handles.dataPcopy = handles.dataP;
    handles.dataPcPcopy = handles.dataPcPorig;
    for i = 1:length(handles.dataPcPcopy)
      t1(i) = getheader(handles.dataPcopy(i),'t1');
      t4(i) = getheader(handles.dataPcPcopy(i),'t4');
      newT1(i) = t1(i) + handles.sumShiftP(i);
      newT4(i) = t4(i) + handles.sumShiftP(i) + handles.sumShiftPcP(i);
      sumPhaseTimeDiff = sumPhaseTimeDiff + (newT4(i) - newT1(i));
    end 
    newTimeAvg = sumPhaseTimeDiff / (length(handles.dataPcPcopy));
    set(handles.text20, 'String',(num2str(newTimeAvg)));

    %% Make map with stations, event, bounce points, and raypaths

    % Use TauP Path to find PcP reflection points at the CMB and the propagation path for each event-station pair; 
    % plotting limits for map are controlled by the event location as well as KNYN latitude and either KNYN or DONT longitude

    clearvars topLatLim bottomLatLim topLonLim bottomLonLim
    % Define latitude limits for map
    topLatLim = (getheader(handles.dataPcP(1),'EVLA')) + 5;
    if ((abs((getheader(handles.dataPcP(1),'EVLO')))) < 60)
      bottomLatLim = -90;
    else 
      bottomLatLim = -76.5;
    end

    % Define longitude limits for map
    if ((getheader(handles.dataPcP(1),'EVLO')) < 0 || (getheader(handles.dataPcP(1),'EVLO')) > 153.33)
      bottomLonLim = 150;
      topLonLim = (getheader(handles.dataPcP(1),'EVLO')) + 20;
    else
      bottomLonLim = abs((getheader(handles.dataPcP(1),'EVLO')) - 20);
      topLonLim = 167;
    end

    % Created map with lat/lon limits as specified above
    axes(handles.axes3);
    axesm('lambertstd','MapParallels',[bottomLatLim topLatLim],'MapLatLimit',[bottomLatLim topLatLim],'MapLonLimit',[bottomLonLim topLonLim]);
    axis off;
    framem on;
    gridm on;
    mlabel on;
    plabel on;
    geoshow('landareas.shp','FaceColor',[0.3 0.3 0.3]);
    hold on;
    guidata(hObject,handles);

    % Read large file of CMB bounce points from the data directory; first two columns are event names and station names - which we 
    % don't need - so they are cleared out immediately; plot these CMB bounce points on map as small grey squares; then clear out big matrices 
    [c1,c2,cmbLats,cmbLons] = textread('./all_evts_CMB_bouncepoints','%s%s%f%f');
    clearvars c1 c2
    geoshow(cmbLats,cmbLons,'DisplayType','point','Marker','s','MarkerSize',3,'MarkerFaceColor',[0.7 0.7 0.7],'MarkerEdgeColor','none');
    clearvars cmbLats cmbLons

    % Use TauP path to determine CMB bounce points and propagation paths; compute average CMB bounce point too
    clearvars bounceLatTot bounceLonTot tt CMB bounceLat bounceLon bounceLatAvg bounceLonAvg baLa baLo BBinfo 
    i = 1;
    bounceLatTot = 0;
    bounceLonTot = 0;

    while (i <= length(handles.dataPcP))
      tt = tauppath('mod','prem','dep',(getheader(handles.dataPcP(i),'EVDP')),'ph','PcP','sta',[(getheader(handles.dataPcP(i),'STLA')),(getheader(handles.dataPcP(i),'STLO'))],'evt',[(getheader(handles.dataPcP(i),'EVLA')),(getheader(handles.dataPcP(i),'EVLO'))]);
      CMB = max(tt.path.depth);
 
      % Plot event-station paths on map
      geoshow(tt.path.latitude(:,1),tt.path.longitude(:,1),'DisplayType','line','Color','k');

      % Find PcP bounce point for this event-station pair and save associated lat/lon
      m = 1;
      while (m <= length(tt.path.depth))
         if (tt.path.depth(m,1) == CMB)
           bounceLat(i) = tt.path.latitude(m,1);
           bounceLon(i) = tt.path.longitude(m,1);
         end
         m = m + 1;
      end

      % Plot CMB bounce points (for this event) and stations on map
      geoshow(bounceLat(i),bounceLon(i),'DisplayType','point','Marker','s','MarkerSize',9,'MarkerFaceColor','c','MarkerEdgeColor','k');
      geoshow((getheader(handles.dataPcP(i),'STLA')),(getheader(handles.dataPcP(i),'STLO')),'DisplayType','point','Marker','^','MarkerSize',9,'MarkerFaceColor','y','MarkerEdgeColor','k');

      % Sum to compute average bounce point location
      bounceLatTot = bounceLatTot + bounceLat(i);
      bounceLonTot = bounceLonTot + bounceLon(i);

      clearvars tt CMB
      i = i + 1;
    end

    % Plot the event on map
    geoshow((getheader(handles.dataPcP(1),'EVLA')),(getheader(handles.dataPcP(1),'EVLO')),'DisplayType','point','Marker','o','MarkerSize',10,'MarkerFaceColor','m','MarkerEdgeColor','k')
    
    % Add average bounce point information to top of GUI window
    bounceLatAvg = bounceLatTot / length(handles.dataPcP);
    bounceLonAvg = bounceLonTot / length(handles.dataPcP);
    baLa = num2str(bounceLatAvg);
    baLo = num2str(bounceLonAvg);
    BBinfo = strcat('Avg CMB Bounce Point  Lat:',32,baLa,32,'Lon:',32,baLo);
    set(handles.text3, 'String', BBinfo);
    guidata(hObject,handles);

    % If want to print to eps file, use following
    outFileeps = strcat(handles.eventname,'.eps');
    set(gcf,'PaperPositionMode','auto');
    saveas(gcf,(sprintf('%s',outFileeps{:})),'epsc');

    % Save the components to a file called evt#.mat; as with GUI, the # refers to the event number in the master list
    % Note that the saved data files (evt.dataP and evt.dataPcP) have all the same header info included (e.g. EVLA,EVO,
    % EVDP, MAG, GCARC,DELTA,KSTNM,T1,T4,etc.) so these don't need to be saved independently.
    evt.event = handles.eventname;
    evt.signalStart = signalStart;
    evt.signalEnd = signalEnd;
    evt.noiseStart = noiseStart;
    evt.noiseEnd = noiseEnd;
    evt.signalStartStack = signalStartStack;
    evt.signalEndStack = signalEndStack;
    evt.noiseStartStack = noiseStartStack;
    evt.noiseEndStack = noiseEndStack;
    evt.searchWindow = limitTimeP/sampRate;
    evt.windowSize = windowSizeT;
    evt.maxShift = maxShift;
    evt.codaAllow = codaAllow;
    evt.stationW = handles.stationW;
    evt.dataP = handles.dataP;
    evt.timeP = timeP;
    evt.sumShiftP = handles.sumShiftP;
    evt.snrP = handles.snrP;
    evt.SNRweightP = SNRweightP;
    evt.coefP = handles.coefP;
    evt.evtWP = handles.evtWP;
    evt.stdDevP = stdDevP;
    evt.weightPstack = handles.weightPstack;
    evt.maxAmpP = [timeMaxP,maxampP];
    evt.PstackSNR = PstackSNR;
    evt.dataPcP = handles.dataPcPorig;
    evt.timePcP = timePcP;
    evt.sumShiftPcP = handles.sumShiftPcP;
    evt.snrPcP = handles.snrPcP;
    evt.snrPcP2P = handles.snrPcP2P;
    evt.SNRweightPcP = SNRweightPcP;
    evt.coefPcP = handles.coefPcP;
    evt.evtWPcP = handles.evtWPcP;
    evt.stdDevPcP = stdDevPcP;
    evt.weightPcPstack = handles.weightPcPstack;
    evt.weightPcPstackPad = handles.weightPcPstackPad;
    evt.maxAmpPcP = [timeMaxPcP,maxampPcP];
    evt.PcPstackSNR = PcPstackSNR;			
    evt.remainder = diffAmp;
    evt.ARstack = ARstack;
    evt.ampRatio = maxamp;
    evt.cmbBPlat = bounceLat;
    evt.cmbBPlon = bounceLon;
    evt.cmbBPlatAvg = bounceLatAvg;
    evt.cmbBPlonAvg = bounceLonAvg;
    evt.PREMavg = PREMavg;
    evt.ShiftTimeAvg = newTimeAvg;
    evt.DiffPcPminusP = (PREMavg - newTimeAvg);   
   
    outFilemat = strcat(handles.eventname,'.mat');
    save((sprintf('%s',outFilemat{:})),'evt');

    % Clear the "saving" array (evt) before proceeding
    clear evt

end

disp('Finished running list of events');

