相对角距离方法的Matlab实现

clear

tic

height=11.1087999344000004;width=11.1087999344000004;long=11.1087999344000004;

load('data.mat')

N=52;

n=60;  %划分区间个数

rc = sqrt(width^2 height^2 long^2); %搜索圆的最大半径

dr = rc/n;%确定半径区间

num = round(rc/dr);

[m,n]=size(a);

%读数据晶型区间划分，每32个坐标数据为一个结构，储存在元胞中

for i=1:length(a)/N

   b= a(1 (N*(i-1)):N*i,:);

   Centers{i,1}=b;

end

%划分距离区间

gr=zeros(num,1);

for n=1:length(Centers)

    centers=Centers{n,1};

    particle_num=length(centers);

    [row,col] = size(centers);

    for i=1:(row-1)

        for j = i 1:row

            %这两个语句的目的是剔除掉原子i,j的坐标

            centersij=centers;

            centersij([i,j],:)=[];

            %增加的判据语句RAD

            a =sqrt(sum((centers(i,:)-centers(j,:)).^2));

            b=sqrt((centers(j,1)-centersij(:,1)).^2 (centers(j,2)-centersij(:,2)).^2 (centers(j,3)-centersij(:,3)).^2);

            c=sqrt((centers(i,1)-centersij(:,1)).^2 (centers(i,2)-centersij(:,2)).^2 (centers(i,3)-centersij(:,3)).^2);

            %b=sqrt(arrayfun(@(a)((centers(j,1)-a).^2),centersij(:,1)) arrayfun(@(a)((centers(j,2)-a).^2),centersij(:,2)) arrayfun(@(a)((centers(j,3)-a).^2),centersij(:,3)));

            %c=sqrt(arrayfun(@(a)((centers(i,1)-a).^2),centersij(:,1)) arrayfun(@(a)((centers(i,2)-a).^2),centersij(:,2)) arrayfun(@(a)((centers(i,3)-a).^2),centersij(:,3)));

            cosa=(b.^2 c.^2-a.^2)./(2*b.*c);

            RAD=1/a.^2-cosa./b.^2;

            distance =sqrt(sum((centers(i,:)-centers(j,:)).^2));

            if distance <= rc%计算

                lane = round(distance/dr);%将粒子划分不同梯度内

                %进行计数前加一个判据条件

                if  RAD>=0 

                    gr(lane) = gr(lane) 1;%做个数累计

                else

                    continue;

                end

            end

        end

    end

end

%绘制径向分布函数

[row,col] = size(gr);

percent = zeros(row,1);

gobal_rho = particle_num / (3/4*rc^3);%全局数密度

for i=1:row

    temp = gr(i);%不同半径下的单位原内的原子个数

    temp = temp / (particle_num*length(Centers));

    temp = temp / (4/3*pi*((i*dr)^3-((i-1)*dr)^3));%某一个半径梯度下的局部密度

    percent(i) = temp / gobal_rho;

end

index=nonzeros(linspace(0,rc,num 1))';

percent = reshape(percent,1,row);

xlim([0, 20]);

ylim([0, 20]);

values = spcrv([[index(1) index index(end)];[percent(1) percent percent(end)]],3)';

r=values(:,1);g=values(:,2);

figure1 = figure;

axes1 = axes('Parent',figure1);

hold(axes1,'on');

plot(index,percent,'LineWidth',3,'Color',[1 0 0]);

plot(r,g,'LineWidth',3,'Color',[0 0 1]);

plot(index,percent,'MarkerFaceColor',[0 1 0],'MarkerSize',8,'Marker','o','LineStyle','none','Color',[0 1 0]);

box(axes1,'on');

hold(axes1,'off');

set(axes1,'FontName','Times New Roman','FontSize',15,'LineWidth',3);

%% 配位数计算

load('data.mat');

fRAD= fit(r,g, 'smoothingspline');

rgRAD=linspace(r(1),r(end),1000)';

intgtotal_cal=rho*4*pi*rgRAD.^2.*fRAD(rgRAD);

z=(max(rgRAD)-min(rgRAD)).*(sum(intgtotal_cal))./1000;