% This script uses the .mat files (for events specified in the list) and assembles all the time shifts
% together.  It then cycles through with each station being the "reference" station and finds the inter-
% station time difference between each station and the reference.  Then it finds the average and standard
% deviation of each station's interstation time differences (again, with each station being the reference).
% And finally it computes the (interstation time difference) - (avg time difference) for each wavefrom
% (for each earthquake).  The results are used to compute corresponding station and event weights
% (seperate event weights for P and PcP), which will then be fed into the MATLAB script to make updated
% stacks with these weights applied.

% Specify input file (list of events)
list = importdata('eventList');
k = 1;

% First, go in and setup a column that lists the station names in alphabetical order; prep TS and TS2
TS(1,1) = {'APRL'};
TS(2,1) = {'BEBP'};
TS(3,1) = {'CASY'};
TS(4,1) = {'DONT'};
TS(5,1) = {'DUBY'};
TS(6,1) = {'FOOT'};
TS(7,1) = {'GRAW'};
TS(8,1) = {'KNYN'};
TS(9,1) = {'LEON'};
TS(10,1) = {'MICH'};
TS(11,1) = {'RAPH'};
TS(12,1) = {'RKST'};
TS(13,1) = {'SAMH'};
TS(14,1) = {'SHRD'};
TS(15,1) = {'SPLN'};
TS2 = TS;

% Go through each .mat file in list and pull the time shifts for the corresponding station.
% If that station doesn't exist for a given event, set time shift value to "NaN"
while (k <= length(list))
  clear event evt
  event = strcat(list(k),'.mat');
  load (sprintf('%s',event{:}));
  q = 1;

  while (q <= size(TS,1))
    found = 0;
    j = 1;
    while (j <= length(evt.sumShiftP))
      if (strcmp((strcat(TS(q,1),'.Pf')),(cellstr(evt.dataP(j).name))))
        % found a matching station entry
        TS(q,k+1) = {evt.sumShiftP(j)};  
        TS2(q,k+1) = {evt.sumShiftPcP(j)};
        found = 1;
      end
      j = j + 1;
    end

    if (found == 0)
      % No matching station; set to "NaN"
      TS(q,k+1) = {NaN};
      TS2(q,k+1) = {NaN};
    end

    q = q + 1;
  end

  k = k + 1;
end

% Now we have all the time shifts assembled into one matrix and organized by station (both for P and PcP).
% A given column (starting at column 2 in TS or TS2) corresponds to the event in the input list.
% We want to cycle through the list of stations, setting each one as the reference station,
% and determine the time shift difference between each station and the corresponding reference 

masterNum = 1;
while (masterNum <= size(TS,1))
  clear timeDiff
  j = 2;
  while (j <= size(TS,2))
    q = 1;
    while (q <= size(TS,1))
      timeDiff(q,j-1) = cell2mat(TS(q,j)) - cell2mat(TS(masterNum,j));
      q = q + 1;
    end
    j = j + 1;
  end

  % Save into structure; number reflects station (1 = APRL, 2 = BEBP, 3 = CASY, etc.)
  masterP(masterNum).timeDiffMS = timeDiff; 

  masterNum = masterNum + 1;
end

% Repeat above for PcP

masterNum = 1;
while (masterNum <= size(TS2,1))
  clear timeDiff
  j = 2;
  while (j <= size(TS2,2))
    q = 1;
    while (q <= size(TS2,1))
      timeDiff(q,j-1) = cell2mat(TS2(q,j)) - cell2mat(TS2(masterNum,j));
      q = q + 1;
    end
    j = j + 1;
  end

  % Save into structure; number reflects station (1 = APRL, 2 = BEBP, 3 = CASY, etc.)
  masterPcP(masterNum).timeDiffMS = timeDiff; 

  masterNum = masterNum + 1;
end

% Now, for each reference set, find the average and standard deviation of the time shifts
% for each station (ignoring any NaN entries)

masterNum = 1;
while (masterNum <= size(TS,1))
  i = 1;
  while (i <= (size(masterP(masterNum).timeDiffMS,1)))
    j = 1;
    count = 0;
    clear pluck;
    r = 1;
    while (j <= (size(masterP(masterNum).timeDiffMS,2)))
      if (isnan(masterP(masterNum).timeDiffMS(i,j)))
        % Value not assigned (station missing); skip
      else
        % extract this value (along this station row)
        pluck(r) = masterP(masterNum).timeDiffMS(i,j);
        r = r + 1;
        count = count + 1;
      end
      j = j + 1;
    end

    if (count == 0)
      % If no events were recorded by this station, set its mean and std dev to NaN
      masterP(masterNum).avgTimeDiff(i,1) = NaN;
      masterP(masterNum).stdDevTimeDiff(i,1) = NaN;
    else
      % Otherwise, compute the mean and std dev for all events at this station, ignoring NaN values
      masterP(masterNum).avgTimeDiff(i,1) = mean(pluck);
      masterP(masterNum).stdDevTimeDiff(i,1) = std(pluck);
    end
    masterP(masterNum).numEvts(i,1) = count;

    i = i + 1;
  end
  masterNum = masterNum + 1;
end

% Repeat above for PcP

masterNum = 1;
while (masterNum <= size(TS2,1))
  i = 1;
  while (i <= (size(masterPcP(masterNum).timeDiffMS,1)))
    j = 1;
    count = 0;
    clear pluck;
    r = 1;
    while (j <= (size(masterPcP(masterNum).timeDiffMS,2)))
      if (isnan(masterPcP(masterNum).timeDiffMS(i,j)))
        % Value not assigned (station missing); skip
      else
        % extract this value (along this station row)
        pluck(r) = masterPcP(masterNum).timeDiffMS(i,j);
        r = r + 1;
        count = count + 1;
      end
      j = j + 1;
    end

    if (count == 0)
      % If no events were recorded by this station, set its mean and std dev to NaN
      masterPcP(masterNum).avgTimeDiff(i,1) = NaN;
      masterPcP(masterNum).stdDevTimeDiff(i,1) = NaN;
    else
      % Otherwise, compute the mean and std dev for all events at this station, ignoring NaN values
      masterPcP(masterNum).avgTimeDiff(i,1) = mean(pluck);
      masterPcP(masterNum).stdDevTimeDiff(i,1) = std(pluck);
    end
    masterPcP(masterNum).numEvts(i,1) = count;

    i = i + 1;
  end
  masterNum = masterNum + 1;
end

% Using the computed mean values, subtract them from the "timeDiffMS" values for each reference set to see how they vary

masterNum = 1;
while (masterNum <= size(TS,1))
  j = 1;
  while (j <= (size(masterP(masterNum).timeDiffMS,2)))
    i = 1;
    while (i <= (size(masterP(masterNum).timeDiffMS,1)))
      masterP(masterNum).timeDiffAvg(i,j) = masterP(masterNum).timeDiffMS(i,j) - masterP(masterNum).avgTimeDiff(i,1);
      i = i + 1;
    end
    j = j + 1;
  end
  masterNum = masterNum + 1;
end

% Repeat above for PcP

masterNum = 1;
while (masterNum <= size(TS2,1))
  j = 1;
  while (j <= (size(masterPcP(masterNum).timeDiffMS,2)))
    i = 1;
    while (i <= (size(masterPcP(masterNum).timeDiffMS,1)))
      masterPcP(masterNum).timeDiffAvg(i,j) = masterPcP(masterNum).timeDiffMS(i,j) - masterPcP(masterNum).avgTimeDiff(i,1);
      i = i + 1;
    end
    j = j + 1;
  end
  masterNum = masterNum + 1;
end

%% Combine values and write out to ascii files %% 

% First, print out all original time shift values (from the catalog plots and mat files)
% (not including station names).
i = 1;
while (i <= size(TS,1))
  j = 2;
  while (j <= size(TS,2))
    TSentriesP(i,j-1) = cell2mat(TS(i,j));
    j = j + 1;
  end
  i = i + 1;
end
dlmwrite('TS_P',TSentriesP,'delimiter','\t');

% Repeat for TS2 (PcP)
i = 1;
while (i <= size(TS2,1))
  j = 2;
  while (j <= size(TS2,2))
    TSentriesPcP(i,j-1) = cell2mat(TS2(i,j));
    j = j + 1;
  end
  i = i + 1;
end
%dlmwrite('TS_PcP',TSentriesPcP,'delimiter','\t');

% Second, cycle through and combine timeDiffMS and timeDiffAvg entries into one big matrix to save to ascii files
masterNum = 1;
k = 1;
while (masterNum <= size(TS,1))
  i = 1;
  while (i <= (size(masterP(masterNum).timeDiffMS,1)))
    j = 1;
    while (j <= (size(TS,2)-1))
      PcombinedTimeDiffMS(k,j) = masterP(masterNum).timeDiffMS(i,j);
      PcombinedTimeDiffAvg(k,j) = masterP(masterNum).timeDiffAvg(i,j);
      j = j + 1;
    end
    k = k + 1;
    i = i + 1;
  end
  masterNum = masterNum + 1;
end
dlmwrite('timeDiffMS_P',PcombinedTimeDiffMS,'delimiter','\t');
dlmwrite('timeDiffAvg_P',PcombinedTimeDiffAvg,'delimiter','\t');

% Repeat above for PcP
masterNum = 1;
k = 1;
while (masterNum <= size(TS2,1))
  i = 1;
  while (i <= (size(masterPcP(masterNum).timeDiffMS,1)))
    j = 1;
    while (j <= (size(TS2,2)-1))
      PcPcombinedTimeDiffMS(k,j) = masterPcP(masterNum).timeDiffMS(i,j);
      PcPcombinedTimeDiffAvg(k,j) = masterPcP(masterNum).timeDiffAvg(i,j);
      j = j + 1;
    end
    k = k + 1;
    i = i + 1;
  end
  masterNum = masterNum + 1;
end
%dlmwrite('timeDiffMS_PcP',PcPcombinedTimeDiffMS,'delimiter','\t');
%dlmwrite('timeDiffAvg_PcP',PcPcombinedTimeDiffAvg,'delimiter','\t');

% Then, combine station-specific average time differences (avgTimeDiff) and std devs (stdDevTimeDiff) into one matrix together.
masterNum = 1;
while (masterNum <= size(TS,1))
  combinedAvgP(:,masterNum) = masterP(masterNum).avgTimeDiff;
  combinedStdDevP(:,masterNum) = masterP(masterNum).stdDevTimeDiff;
  masterNum = masterNum + 1;
end
dlmwrite('timeShiftAvgs_P',combinedAvgP,'delimiter','\t');
dlmwrite('timeShiftStdDev_P',combinedStdDevP,'delimiter','\t');

% Repeat above for PcP
masterNum = 1;
while (masterNum <= size(TS2,1))
  combinedAvgPcP(:,masterNum) = masterPcP(masterNum).avgTimeDiff;
  combinedStdDevPcP(:,masterNum) = masterPcP(masterNum).stdDevTimeDiff;
  masterNum = masterNum + 1;
end
%dlmwrite('timeShiftAvgs_PcP',combinedAvgPcP,'delimiter','\t');
%dlmwrite('timeShiftStdDev_PcP',combinedStdDevPcP,'delimiter','\t');

%% Compute weights %%

%% FIRST - STATION WEIGHTS %%

% Using the P-wave standard deviations (i.e., combinedStdDevP), compute the associated
% L2 norm value for each column (each station as reference)

j = 1;
while (j <= (size(combinedStdDevP,2)))
  sqsum = 0;
  i = 1;
  while (i <= (size(combinedStdDevP,1)))
    sqsum = sqsum + (combinedStdDevP(i,j)^2);
    i = i + 1;
  end
  L2norm(j) = sqsum^0.5;
  j = j + 1;
end

% Find min and max L2norm value
minL2norm = min(L2norm);
maxL2norm = max(L2norm);

% Compute fraction of each L2norm value compared to min and max values
% and use that to compute a station weight (specifying the maximum amount
% of downweighting to use for the "worst" station)
i = 1;
maxDownWeight = 30;

while (i <= length(L2norm))
  L2normFract(i) = (L2norm(i) - minL2norm) / (maxL2norm - minL2norm);
  stationWeight(i) = (100 - (maxDownWeight*L2normFract(i)))/ 100;
  i = i + 1;
end
dlmwrite('stationWeights',stationWeight','delimiter','\t');


%% NEXT - EVENT WEIGHTS %%

% Using the P-wave differences from the averages (i.e., PcombinedTimeDiffAvg), compute the
% L2 norm values for each station-event 
j = 1;
while (j <= size(PcombinedTimeDiffAvg,2))
  i = 1;
  k = 1;
  sqsum = 0;
  while (i <= size(PcombinedTimeDiffAvg,1))
    if (i == 16 || i == 31 || i == 46 || i == 61 || i == 76 || i == 91 || i == 106 || i == 121 || i == 136 || i == 151 || i == 166 || i == 181 || i == 196 || i == 211) 
      L2normP(k,j) = sqsum ^ 0.5;
      k = k + 1;
      sqsum = 0;
    end

    if (isnan(PcombinedTimeDiffAvg(i,j)) == 0)
      sqsum = sqsum + (PcombinedTimeDiffAvg(i,j)^2);
    end

    if (i == 225)
      L2normP(k,j) = sqsum ^ 0.5;
    end

    i = i + 1;
  end
  j = j + 1;
end

% Compute mean L2 norm value and (0.5 * std dev) L2 norm value for each station (i.e., each row in L2norm matrix created above)
i = 1;
while (i <= size(L2normP,1))
  L2normAvgP(i) = mean(L2normP(i,:));
  L2normStdDevP(i) = (std(L2normP(i,:)))*0.5;
  i = i + 1;
end

% Write out average and standard deviation L2 norm values for later use
dlmwrite('L2normAvg_P',L2normAvgP,'delimiter','\t');
dlmwrite('L2normStdDev_P',L2normStdDevP,'delimiter','\t');

% Compute associated weight
i = 1;
while (i <= size(L2normP,1))
  j = 1;
  while (j <= size(L2normP,2))
    if (L2normP(i,j) == 0)
      eventWeightP(i,j) = NaN;
    elseif (L2normP(i,j) < (L2normAvgP(i) + L2normStdDevP(i)))
      eventWeightP(i,j) = 1;
    else
      eventWeightP(i,j) = 1 - ((L2normP(i,j)-(L2normAvgP(i)+L2normStdDevP(i))) / L2normAvgP(i));
      if (eventWeightP(i,j) < 0)
        eventWeightP(i,j) = 0;
      end
    end
    j = j + 1;
  end
  i = i + 1;
end
dlmwrite('eventWeights_P',eventWeightP,'delimiter','\t');

% Report how many event-station records are getting a weight less than 1 across all data (not including NaN values)
% and reform into vector (without NaN values) to make a histogram plot
i = 1;
totcount = 0;
weightcount = 0;
while (i <= size(eventWeightP,1))
  j = 1;
  while (j <= size(eventWeightP,2))
    if (isnan(eventWeightP(i,j)) == 0)
      totcount = totcount + 1;
      H1(totcount) = eventWeightP(i,j);
      if (eventWeightP(i,j) < 1)
        weightcount = weightcount + 1;
        H2(weightcount) = eventWeightP(i,j);
      end
    end
    j = j + 1;
  end
  i = i + 1;
end
message1 = strcat('Total number of P-wave event-station weights:',32,num2str(totcount));
message2 = strcat('Total number of P-wave weights less than one:',32,num2str(weightcount));
disp(message1);
disp(message2);

figure(1);
subplot(2,1,1);
hist(H1,10);
title('P-waves, All Weights')
xlabel('weight')
ylabel('count')
subplot(2,1,2);
hist(H2);
title('P-waves, Weights < 1')
xlabel('weight')
ylabel('count')

% Repeat above for PcP event weights
j = 1;
while (j <= size(PcPcombinedTimeDiffAvg,2))
  i = 1;
  k = 1;
  sqsum = 0;
  while (i <= size(PcPcombinedTimeDiffAvg,1))
    if (i == 16 || i == 31 || i == 46 || i == 61 || i == 76 || i == 91 || i == 106 || i == 121 || i == 136 || i == 151 || i == 166 || i == 181 || i == 196 || i == 211) 
      L2normPcP(k,j) = sqsum ^ 0.5;
      k = k + 1;
      sqsum = 0;
    end

    if (isnan(PcPcombinedTimeDiffAvg(i,j)) == 0)
      sqsum = sqsum + (PcPcombinedTimeDiffAvg(i,j)^2);
    end

    if (i == 225)
      L2normPcP(k,j) = sqsum ^ 0.5;
    end

    i = i + 1;
  end
  j = j + 1;
end

% Compute mean L2 norm value and (0.5 * std dev) L2 norm value for each station (i.e., each row in L2norm matrix created above)
i = 1;
while (i <= size(L2normPcP,1))
  L2normAvgPcP(i) = mean(L2normPcP(i,:));
  L2normStdDevPcP(i) = (std(L2normPcP(i,:)))*0.5;
  i = i + 1;
end

% Write out average and standard deviation L2 norm values for later use
%dlmwrite('L2normAvg_PcP',L2normAvgPcP,'delimiter','\t');
%dlmwrite('L2normStdDev_PcP',L2normStdDevPcP,'delimiter','\t');

% Compute associated weight
i = 1;
while (i <= size(L2normPcP,1))
  j = 1;
  while (j <= size(L2normPcP,2))
    if (L2normPcP(i,j) == 0)
      eventWeightPcP(i,j) = NaN;
    elseif (L2normPcP(i,j) < (L2normAvgPcP(i) + L2normStdDevPcP(i)))
      eventWeightPcP(i,j) = 1;
    else
      eventWeightPcP(i,j) = 1 - ((L2normPcP(i,j)-(L2normAvgPcP(i)+L2normStdDevPcP(i))) / L2normAvgPcP(i));
      if (eventWeightPcP(i,j) < 0)
        eventWeightPcP(i,j) = 0;
      end
    end
    j = j + 1;
  end
  i = i + 1;
end
%dlmwrite('eventWeights_PcP',eventWeightPcP,'delimiter','\t');

% Report how many event-station records are getting a weight less than 1 across all data (not including NaN values)
% and reform into vector (without NaN values) to make a histogram plot
i = 1;
totcount = 0;
weightcount = 0;
while (i <= size(eventWeightPcP,1))
  j = 1;
  while (j <= size(eventWeightPcP,2))
    if (isnan(eventWeightPcP(i,j)) == 0)
      totcount = totcount + 1;
      H3(totcount) = eventWeightPcP(i,j);
      if (eventWeightPcP(i,j) < 1)
        weightcount = weightcount + 1;
        H4(weightcount) = eventWeightPcP(i,j);
      end
    end
    j = j + 1;
  end
  i = i + 1;
end
message3 = strcat('Total number of PcP-wave event-station weights:',32,num2str(totcount));
message4 = strcat('Total number of PcP-wave weights less than one:',32,num2str(weightcount));
disp(message3);
disp(message4);

figure(2);
subplot(2,1,1);
hist(H3,10);
title('PcP-waves, All Weights')
xlabel('weight')
ylabel('count')
subplot(2,1,2);
hist(H4);
title('PcP-waves, Weights < 1')
xlabel('weight')
ylabel('count')


