%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: BlockBased3D_RangeSearch(puctr,puradius,dsites,index)
%
% Goal: find the data sites located in a given subdomain and the distances
%       between the subdomain centre and data sites
%
% Inputs: puctr:     subdomain centre
%         puradius:  radius of PU subdomains
%         dsites:    NX3 matrix representing a set of N data sites
%         index:     vector containing the indices of the data points 
%                    located in the k-th block (the cell containing the 
%                    subdomain centre) and in the neighbouring cells
%  
% Outputs: idx:  vector containing the indices of the data points located
%                in a given PU subdomain
%          dist: vector containing the distances between the data sites 
%                and the subdomain centre
%
%-------------------------------------------------------------------------%
function [idx,dist] = BlockBased3D_RangeSearch(puctr,puradius,dsites,index)
N = size(dsites,1); dist1 = zeros(1,N); dist = []; idx = []; %Initialize
% Compute distances between the data sites and the centre
for i = 1:N
    dist1(i) = sqrt((puctr(1,1)-dsites(i,1))^2 + ...
        (puctr(1,2)-dsites(i,2))^2 + (puctr(1,3)-dsites(i,3))^2);
end
% Use a sort procedure to order the distances
[sort_dist,IX] = sort(dist1);
N1 = size(sort_dist,2); j1 = 1; j2 = 1; %Initialize
% Find the data sites located in a given subdomain
while (j2 <= N1) && (sort_dist(j2) <= puradius)
        if nargin == 3
            idx(j1) = IX(j2); dist(j1) = dist1(IX(j2)); 
            j1 = j1 + 1; j2 = j2 + 1;
        else
            idx(j1) = index(IX(j2)); dist(j1) = dist1(IX(j2));
            j1 = j1 + 1; j2 = j2 + 1;
        end
end