function [ref_match, twin_match] = gt6DGetMatchingReflections(ref_or, ors, angular_toll, cryst)
    % ref_match contains the indices of shared reflections in the
    % twin_gr.allblobs list: i.e. [0 10; ...] would tell us that
    % ref.allblobs(1) shares a spot with twin_or(2).allblobs(10)
    %
    % twin_match contains the indices of shared reflections in the
    % ref_gr.allblobs list: i.e. [10 0; ...] would tell us that
    % twin_or(1).allblobs(1) shares a spot with ref_gr.allblobs(10)

    if (~exist('angular_toll', 'var') || isempty(angular_toll))
        angular_toll = 1; % degrees
    end

    ref_ws = ref_or.allblobs.omega;
    ref_ns = ref_or.allblobs.eta;
    ref_tt = ref_or.allblobs.thetatype;
    ref_Lfactors = ref_or.allblobs.lorentzfac;

    w_toll = angular_toll * ref_Lfactors;
    n_toll = angular_toll;

    ors_allblobs = [ors(:).allblobs];
    ors_ws = cat(2, ors_allblobs(:).omega);
    ors_ns = cat(2, ors_allblobs(:).eta);
    ors_tt = cat(2, ors_allblobs(:).thetatype);

    num_ors = numel(ors);
    num_spots = numel(ref_ws);
    ones_array = ones(num_spots, num_ors);
    ref_match = zeros(num_spots, num_ors);
    twin_match = zeros(num_spots, num_ors);

    row_ind = [];

    for ii = 1 : num_spots
        ref_ws_array = ref_ws(ii) * ones_array;
        ref_ns_array = ref_ns(ii) * ones_array;
        ref_tt_array = ref_tt(ii) * ones_array;
        tt_equal = ref_tt_array == ors_tt;
        w_diff = abs(mod(ref_ws_array - ors_ws + 180, 360) - 180);
        n_diff = abs(mod(ref_ns_array - ors_ns + 180, 360) - 180);
        [row_ind, col_ind] = find((w_diff < w_toll(:,ones(num_ors,1))) & (n_diff < n_toll) & tt_equal);
        if ~isempty(row_ind)
            for jj = 1 : numel(row_ind)
                ref_match(ii, col_ind(jj)) = row_ind(jj);
                twin_match(row_ind(jj), col_ind(jj)) = ii;
            end
        end
    end

end