[xyzd, Axons_Number, Bin_Size]=axon_dist(Path_FileName,yz)

function [xyzd, Axons_Number, Bin_Size]=axon_dist(Path_FileName,yz)
    FID=fopen(Path_FileName,'r');
    line=['zeng'];
    Axons_Distribution=[];
    while isempty(str2num(line))
        line=fgetl(FID);
    end
    while line~=-1
        Axons_Distribution=[Axons_Distribution;str2num(line)];
        line=fgetl(FID);
    end
    fclose(FID);
    if size(Axons_Distribution,1)==1
        Bin_Size=[];
    else
        temp=diff(Axons_Distribution(:,1));
        if ~isempty(find(temp(2:length(temp))~=temp(1)))
            disp('Bin size is not uniform')
        end
        Bin_Size=temp(1);
    end
    Axons_Number=sum(Axons_Distribution(:,2));
    locations=length(yz(:,1));
    if Axons_Number>locations
   	    zexst('err','this distribution has too many axons')   
        xyzd=0;
        Axons_Number=0;
    else
        axon_orders=randperm(locations);
        axon_orders=axon_orders(1:Axons_Number)';
        Dia_Distribution=[];
        for i = 1:length(Axons_Distribution(:,1))
            Dia_Distribution=[Dia_Distribution;Axons_Distribution(i,1)*ones(Axons_Distribution(i,2),1)];
        end
        %For histrogram at the center of bins. It is a must; other,
        %Plot_Range at "Apply" will be a problem
        Dia_Distribution=Dia_Distribution+Bin_Size*(round(rand(Axons_Number,1)*99)/100-0.5);
        xyzd=zeros(Axons_Number,4);
        xyzd(:,[2 3])=yz(axon_orders,:);
        xyzd(:,[1 4])=[rand([Axons_Number 1])-0.5 Dia_Distribution];
    end
   
    %------------------------------------
    if 0
        y=600;
        z=200;
%        z=200;
        temp=find(xyzd(:,2)==y & xyzd(:,3)==z);
        if isempty(temp)
            disp(['XYZD(1) is manually assigned in Axon_Dist'])
            xyzd(1,:)=[0 y z 5];
%            xyzd(1,:)=[1 y z 7.5];
%            xyzd(1,:)=[0.5 y z 15];
%            xyzd(1,:)=[0.125 y z 20];
    end
end