x=xlsread('jiagong.xls','I2:I282');
y=xlsread('jiagong.xls','J2:J282');

disxy=zeros(281,281);
disxyPX=zeros(281,281);
for i=1:281
    for j=1:281
        disxy(i,j)=sqrt(((x(i)-x(j))*(111*cos((y(i)+y(j))/2)))^2+((y(i)-y(j))*111)^2);%单位km，坐标WGS1984
    end
end
%排序
for i=1:281
    a=[];
    for j=1:281
        a(j)=disxy(i,j);
    end
    sort(a);
    for j=1:281
        disxy(i,j)=a(j);
    end
end
     
t1=1;
K=zeros(1,20);
K0=zeros(1,20);
L=zeros(1,20);
for t=5:5:100
    % t是距离，单位km
    A=19400.970606; %研究区面积，单位km2
    n0=281; % 点个数
    count=0;
    for m=1:281
        for n=1:281
            if disxy(m,n)<=t
                count=count+1; %距离小于半径的点数
            end
        end
    end
    K(t1)=(A/n0^2)*count;  %A研究区面积，n为研究区内点事件数量
    K0(t1)=pi*(t^2);     %K0<K，集聚分布
    L(t1)=sqrt(K(t1)/pi)-t;   %L>0，集聚分布
    t1=t1+1;
end