function Y=ftdotp0(X,R,a)
i=1;
j=1;
m=1;
N=size(X,1);
r=zeros(1,N);
while (i<N)
    A=(X(i+1,3)-X(i,3))/(X(i+1,1)-X(i,1));
    while (j<N-i-1)
        K=zeros(1,N-i-1);
        K(j)=(X(i+j+1,3)-X(i,3))/(X(j+i+1,1)-X(i,1));
        if abs((K(j)-A)/(1-K(j)*A))<tan(a*pi/180)
            j=j+1;
        else
            r(m)=sqrt(X(j+i,3)-X(i,3).^2+(X(j+i,1)-X(i,1)).^2)
            m=m+1;
            i=i+j;
            j=1;
            break
        end
    end
a=length(find(r>0));
s=0:0.02:R;
Y=zeros(length(s),2);
Y(:,1)=s';
   for i=1:1:length(s)
       Y(i,2)=length(find(r>s(i)))/a;
   end
end


求指点，谢谢