function [rand_ltimes, rand_ptimes, rand_ptimesnum, rand_llatlon, rand_lEnergy, rand_vecl, rand_vec] = rand_peak_times(winltimes, winllatlon, winlEnergy, winpeaksrmstime, dt, peakwint, B_output, siteData, wintr3, wintc3)
global fs
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)
%         % 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) = [];
%     
%     [wintr3,ia,ic] = unique(wintr2);  % 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.
%     wintc3 = wintc2(ia);
    
    % Create vector of random indices same length as matched selected
    % indices length (wintc3)
    rand_vec{i} = sort(randperm(length(winpeaksrmstime{i}{1}),length(wintr3{i}))); % wintc3
    rand_vecl{i} = sort(randperm(length(winltimes{i}),length(wintr3{i}))); % wintr3
    % Re-associate vectors
    rand_ltimes{i} = winltimes{i}(rand_vecl{i});
    
    % Find lightning distances from station based on lat and lon
    STN_Lat=siteData.geoLat; STN_Long=siteData.geoLon-360;
    
    
    % 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_vecl{i},:);
    rand_lEnergy{i} = winlEnergy{i}(rand_vecl{i},:);
    %rand_matched_dt{i} = diag(wintdiffs(rand_vecl{i},rand_vec{i}));
    [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{i})';
    rand_ptimes{i,2} = [];
    rand_pnegtimes = 0;
    
    for j = 1:length(rand_ptimes{i,1})
        if round(rand_ptimes{i,1}(j)-0.4*peakwint*fs) < 1
            %             ptimes{i,3}(j,:) = 0; % spike window
            continue
        else
            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
        rand_ptimes{i,4}(j,:) = rand_ptimes{i,3}(j,:)/abs(B_output(rand_ptimes{i,1}(j)));
        
        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 = rand_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

    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

    rand_ptimesnum(1,i) = length(rand_ptimes{i});
    rand_ptimesnum(2,i) = rand_ptimesnum(1,i) - rand_pnegtimes;
    rand_ptimesnum(3,i) = rand_pnegtimes;
end
end