function horizon=knn3D(coord,delta,dx)
% 领域识别子函数3D
% coord为坐标，delta是近场范围尺寸，dx是网格划分尺寸
% horizon为领域编号矩阵，第一列为域内点个数的总数，后面几列为相应的近场域内点的编号
total_node=size(coord,1); %物质点总数
horizon_number=(2*floor(delta/dx))^3; %预计的近场域内点的个数，假设为周围一个方块
horizon=zeros(total_node,horizon_number); %领域编号矩阵初始化
for i=1:total_node %对所有点做循环
    %先判断是否在矩形区域
    kk=zeros(total_node,1); %初步符合的点的编号集合
    nnn=0; %计数
    for nn=1:total_node  %再对所有点做一次循环
        %点不能重合，各个方向绝对值要小于delta
        if i~=nn&&abs(coord(i,1)-coord(nn,1))<delta&&abs(coord(i,2)-coord(nn,2))<delta&&abs(coord(i,3)-coord(nn,3))<delta
            nnn=nnn+1;
            kk(nnn)=nn;
        end
    end
    jud2= kk~=0; %找到不为零的编号
    k=kk(jud2); %重编
    %再细化为一个圆
    num=1; %计数
    for j=1:1:numel(k)    %对所有的符合上面一个判断的点进行循环
        cnn=k(j); %找到该点编号
        %算两点间距离
        l=sqrt((coord(cnn,1)-coord(i,1)).^2+(coord(cnn,2)-coord(i,2)).^2+(coord(cnn,3)-coord(i,3)).^2);        
        if l<delta
            num=num+1;
            horizon(i,num)=cnn;
        end
    end
    horizon(i,1)=num-1;
end
%对horizon进行修正，略掉后面为零的编号
mm=max(horizon(:,1));
horizon=horizon(:,1:mm+1);
end
