function [wintlims, ltimes, ptimes, ptimesnum, llatlon, ldist, lEnergy, matched_dt, wintr3, wintc3] = matchedtimes(winltimes, winllatlon, winlEnergy, winpeaksrmstime, dt, peakwint, time, B_output, siteData, heightfact2)
global fs
wintlims = [-0.2, 0.8];
winoverlap = 0.1;
for i = 1:length(winltimes)
% Calculate time in seconds of each lightning strike and B peak
% above the threashold
A = 3600*hour(winltimes{i}) + 60*minute(winltimes{i}) + second(winltimes{i}); % N x 1 Column
B = 3600*hour(winpeaksrmstime{i}{2}) + 60*minute(winpeaksrmstime{i}{2}) + second(winpeaksrmstime{i}{2}); % 1 x M row
% Subtract all B peak times from all lightning times to create
% matrix. Columns represent one B peak time subtracted from each
% lightning strike time.
wintdiffs = A - B;
% Logically, lightning must precede a B peak. Therefore only
% negative values in each column associated to one B peak are
% logically causal. Furthermore, a minimum separation time between
% strike and peak has been set as an input, in this function denoted as
% "dt", in higher level functions as "strikepeakdt". Therefore, only
% negative values with an absolute value less than "dt" must be
% considered, based on the input from the user.
% The next part is where it gets tricky. Earlier iterations of this
% "matchedtimes.m" function selected the minimum value that met the
% logical causal criteria (lightning before B peak) and the minimum
% "dt" input by the user. The argument was that it should be the
% strike nearest in time to the B peak that can be correlated.
% However, without considering the distance between that strike and the
% location of the ICM station, coincidence can cause that minimum
% separation time to be far too low for distant strikes to physically
% reach the ICM detector. In other words, earlier iterations of this
% code allowed for information to travel faster than the speed of
% light. Therefore, distance must be considered next after the causal
% and input "dt" criteria. Furthermore, the number of strikes that
% meet the first two criteria, and can physically reach the detector in
% time could still be greater than one. The question then is which one
% correlates?
% First sort by causal and "dt"
%[wintr, wintc] = find(abs(wintdiffs) == min(abs(wintdiffs)) & wintdiffs < 0 & abs(wintdiffs) < dt); % This is BAAAD. By having abs(wintdiffs) == min(abs(wintdiffs)), a NON-NEGATIVE (non causal) case could be the minimum value. That would mean that the causal case, which passes the other two tests, would still FAIL, because it is not the minimum. This reduces the statistics!!! Eliminates perfectly legitimate cases!!!
% The line above determines the row index "wintr" of the smallest
% value that matches all the criteria in each column. "wintc" gives the
% index of all the columns that have a value that matches all three
% criteria. Both are one dimensional vectors.
[wintr1, wintc1] = find(wintdiffs < 0 & abs(wintdiffs) > dt(1) & abs(wintdiffs) < dt(2)); % Determine the...
if isempty(wintr1) && isempty(wintc1)
wintr4{i} = [];
wintc4{i} = [];
ltimes{i} = [];
llatlon{i} = [];
ldist{i} = [];
lEnergy{i} = [];
matched_dt{i} = [];
llatlon{i} = [NaN NaN 0 1];
ptimes{i,1} = []; ptimes{i,2} = []; ptimes{i,3} = []; ptimes{i,4} = []; ptimes{i,5} = []; ptimes{i,6} = []; ptimes{i,7} = []; ptimes{i,8} = [];
ptimesnum(1,i)=0; ptimesnum(2,i)=0; ptimesnum(3,i)=0;
% random versions
% rand_ltimes{i} = [];
% rand_llatlon{i} = [];
% rand_lEnergy{i} = [];
% rand_matched_dt{i} = [];
% rand_llatlon{i} = [NaN NaN 0 1];
% rand_ptimes{i,1} = []; ptimes{i,2} = []; ptimes{i,3} = []; ptimes{i,4} = []; ptimes{i,5} = []; ptimes{i,6} = []; ptimes{i,7} = []; ptimes{i,8} = [];
% rand_ptimesnum(1,i)=0; ptimesnum(2,i)=0; ptimesnum(3,i)=0;
continue
end
% Here need to do additional filtering to select one B for each
% lightning strike, but the actual "dt" between those two must be
% greater than the time it takes for the strike to travel the distance
% to the ICM detector at the speed of light. If there are still
% multiple, this code picks the lowest "dt", ie. the first causal B
% event after the lightning strike. The justification for doing so, as
% well as the nature of the follow on B events that also match that
% strike, is still being discussed (as of June 2018)
% Find lightning distances from station based on lat and lon
STN_Lat=siteData.geoLat; STN_Long=siteData.geoLon-360;
[arclen,az] = distance(STN_Lat,STN_Long,winllatlon{i}(wintr1,1),winllatlon{i}(wintr1,2));
% Calculate minimum travel times for lightning to reach ICM detector
min_l_travel_time = deg2km(arclen)/300000; % in seconds. distances and speed of light in km and km/s
l_peak_dt = second(winpeaksrmstime{i}{2}(wintc1)') - second(winltimes{i}(wintr1));
%l_peak_dt(l_peak_dt < 0.2) = 0;
l_p_c_diff = l_peak_dt - min_l_travel_time;
wintr2 = wintr1;
wintc2 = wintc1;
wintr2(l_p_c_diff(:,1) < 0) = [];
wintc2(l_p_c_diff(:,1) < 0) = [];
[~,i1] = unique(wintr2,'first'); % This line gives the "flat" indices "ia" of all the row indices in "wintr". The purpose is to ensure one lightning strike can only be matched to one B peak.
c = wintc2(i1);
[b,i2] = unique(wintr2,'last'); % By doing this double run, any cases where there is more than one lightning strike matching to one B or vice versa is GONE
wintr3 = b(i1 == i2);
wintc3 = c(i1 == i2);
[~,i3] = unique(wintc3,'first'); % This line gives the "flat" indices "ia" of all the row indices in "wintr". The purpose is to ensure one lightning strike can only be matched to one B peak.
e = wintr3(i3);
[d,i4] = unique(wintc3,'last'); % By doing this double run, any cases where there is more than one lightning strike matching to one B or vice versa is GONE
wintc4{i} = d(i3 == i4);
wintr4{i} = e(i3 == i4);
% Create vector of random indices same length as matched selected
% indices length (wintc3)
% rand_vec = sort(randperm(length(winpeaksrmstime{i}{1}),length(wintc3)));
% Re-associate vectors
ltimes{i} = winltimes{i}(wintr4{i});
% rand_ltimes{i} = winltimes{i}(rand_vec);
% Find lightning distances from station based on lat and lon
STN_Lat=siteData.geoLat; STN_Long=siteData.geoLon-360;
llatlon{i} = winllatlon{i}(wintr4{i},:);
lEnergy{i,1} = winlEnergy{i}(wintr4{i},:);
matched_dt{i,1} = diag(wintdiffs(wintr4{i},wintc4{i}));
[arclen,az] = distance(STN_Lat,STN_Long,llatlon{i}(:,1),llatlon{i}(:,2));
llatlon{i,1} = [llatlon{i} deg2km(arclen) ones(length(ltimes{i}),1)];
ldist{i,1} = deg2km(arclen);
ptimes{i,1} = winpeaksrmstime{i}{1}(wintc4{i})';
ptimes{i,2} = winpeaksrmstime{i}{2}(wintc4{i})';
pnegtimes = 0;
% Create random selection of B peaks, not selected for proximity to
% lightning. Take N same number as the lightning selected number.
% rand_llatlon{i} = winllatlon{i}(rand_vec,:);
% rand_lEnergy{i} = winlEnergy{i}(rand_vec,:);
% rand_matched_dt{i} = diag(wintdiffs(rand_vec,rand_vec));
% [rand_arclen,rand_az] = distance(STN_Lat,STN_Long,rand_llatlon{i}(:,1),rand_llatlon{i}(:,2));
% rand_llatlon{i} = [rand_llatlon{i} deg2km(rand_arclen) ones(length(rand_ltimes{i}),1)];
% rand_ptimes{i,1} = winpeaksrmstime{i}{1}(rand_vec)';
% rand_ptimes{i,2} = [];
% rand_pnegtimes = 0;
k = 1;
win1offset = 0;
for j = 1:length(ptimes{i,1})
if round(ptimes{i,1}(j)+wintlims(1)*peakwint*fs) < 1
% ptimes{i,3}(j,:) = 0; % spike window
win1offset = win1offset + 1;
continue
else
ptimes{i,12}(j-win1offset,:) = time(round(ptimes{i,1}(j)+wintlims(1)*peakwint*fs) : round(ptimes{i,1}(j)+wintlims(2)*peakwint*fs-1));
ptimes{i,3}(j-win1offset,:) = B_output(round(ptimes{i,1}(j)+wintlims(1)*peakwint*fs) : round(ptimes{i,1}(j)+wintlims(2)*peakwint*fs-1)); % spike window
%rand_ptimes{i,3}(j,:) = B_output(round(rand_ptimes{i,1}(j)-0.4*peakwint*fs) : round(rand_ptimes{i,1}(j)+0.6*peakwint*fs-1));
end
% October 10, 2018 adding code to reject windows that have one or
% more valid peak/strike matches within the superposition. Don't
% care about what comes before, so only reject windows that have
% multiples after.
if j>1 + win1offset
if ptimes{i,2}(j) > ptimes{i,12}(j-win1offset-1,end-winoverlap*peakwint*fs)
ptimes{i,13}(k) = j-1;
k = k+1;
end
end
ptimes{i,7}(j-win1offset,:) = B_output(round(ptimes{i,1}(j)+0.6*peakwint*fs) : round(ptimes{i,1}(j)+1.6*peakwint*fs-1));
[M,N] = size(ptimes{i,3});
ptimes{i,4}(j-win1offset,:) = ptimes{i,3}(j-win1offset,:)/abs(B_output(ptimes{i,1}(j))); % Normalization of spike window by the spike in that window
%rand_ptimes{i,4}(j,:) = rand_ptimes{i,3}(j,:)/abs(B_output(rand_ptimes{i,1}(j)));
if B_output(ptimes{i,1}(j)) < 0
ptimes{i,5}(j-win1offset,:) = -ptimes{i,3}(j-win1offset,:); % reversed polarity
ptimes{i,8}(j-win1offset,:) = -ptimes{i,7}(j-win1offset,:);
ptimes{i,6}(j-win1offset,:) = ptimes{i,5}(j-win1offset,:)/abs(B_output(ptimes{i,1}(j))); % Normalization of spike window by the spike in that window
llatlon{i}(j-win1offset,4) = -1;
pnegtimes = pnegtimes + 1;
else
ptimes{i,5}(j-win1offset,:) = ptimes{i,3}(j-win1offset,:);
ptimes{i,6}(j-win1offset,:) = ptimes{i,3}(j-win1offset,:)/abs(B_output(ptimes{i,1}(j)));
end
%if B_output(rand_ptimes{i,1}(j)) < 0
%rand_ptimes{i,5}(j,:) = -rand_ptimes{i,3}(j,:); % reversed polarity
%rand_ptimes{i,8}(j,:) = -rand_ptimes{i,7}(j,:);
% rand_ptimes{i,6}(j,:) = rand_ptimes{i,5}(j,:)/abs(B_output(rand_ptimes{i,1}(j))); % Normalization of spike window by the spike in that window
% rand_llatlon{i}(j,4) = -1;
% rand_pnegtimes = pnegtimes + 1;
%else
% rand_ptimes{i,5}(j,:) = rand_ptimes{i,3}(j,:);
% rand_ptimes{i,6}(j,:) = rand_ptimes{i,3}(j,:)/abs(B_output(rand_ptimes{i,1}(j)));
%end
end
% Create vector of matched_dt for unique cases
matched_dt{i,2} = matched_dt{i,1}(ptimes{i,13});
llatlon{i,2} = llatlon{i,1}(ptimes{i,13});
ldist{i,2} = ldist{i,1}(ptimes{i,13});
lEnergy{i,2} = lEnergy{i,1}(ptimes{i,13});
% Calculate indices of quantiles of lEnergy, create ptimes subset of
% largest ones
% minenergy = quantile(lEnergy{i},heightfact2);
% minenergy_ind = find(lEnergy{i} >= minenergy);
% minenergy_ind_unique = find(lEnergy{i}(ptimes{i,13}) >= minenergy);
ptimes{i,9} = [];%ptimes{i,3}(minenergy_ind-win1offset,:);
ptimes{i,10} = [];%ptimes{i,5}(minenergy_ind-win1offset,:);
ptimes{i,11} = [];%ptimes{i,2}(minenergy_ind);
ptimes{i,14} = [];%ptimes{i,13}(minenergy_ind_unique)-win1offset;
if length(ptimes{i,3}) < peakwint*fs
ptimes{i,3} = []; % Maybe check this condition
end
if length(ptimes{i,4}) < peakwint*fs
ptimes{i,4} = []; % Maybe check this condition
end
if length(ptimes{i,5}) < peakwint*fs
ptimes{i,5} = []; % Maybe check this condition
end
% if length(rand_ptimes{i,3}) < peakwint*fs
% rand_ptimes{i,3} = []; % Maybe check this condition
% end
% if length(rand_ptimes{i,4}) < peakwint*fs
% rand_ptimes{i,4} = []; % Maybe check this condition
% end
% if length(rand_ptimes{i,5}) < peakwint*fs
% rand_ptimes{i,5} = []; % Maybe check this condition
% end
% if M < 2
% ptimes{i,5} = ptimes{i,3};
% ptimes{i,6} = ptimes{i,4};
% ptimes{i,9} = ptimes{i,8};
% else
% ptimes{i,5} = mean(ptimes{i,3});
% ptimes{i,6} = mean(ptimes{i,4});
% ptimes{i,9} = mean(ptimes{i,8});
% end
% if isnan(ptimes{i,5})
% ptimes{i,5} = [];
% end
% if isnan(ptimes{i,6})
% ptimes{i,6} = [];
% end
% if isnan(ptimes{i,9})
% ptimes{i,9} = [];
% end
ptimesnum(1,i) = length(ptimes{i});
ptimesnum(2,i) = ptimesnum(1,i) - pnegtimes;
ptimesnum(3,i) = pnegtimes;
ptimesnum(4,i) = length(ptimes{i,13});
ptimesnum(5,i) = ptimesnum(4,i) - length(find(B_output(ptimes{i,1}(ptimes{i,13})) < 0));
ptimesnum(6,i) = length(find(B_output(ptimes{i,1}(ptimes{i,13})) < 0));
% rand_ptimesnum(1,i) = length(ptimes{i});
% rand_ptimesnum(2,i) = ptimesnum(1,i) - pnegtimes;
% rand_ptimesnum(3,i) = pnegtimes;
end
end