求球面均布spiral point的matlab程序

我按照别的论文里面的算法，用matlab试了好多次都不合适，麻烦各位了，

你好，这个算法很多人都实现过了，我之前也用静电场作用方法来获取球面均匀分布的点。可以参考一下：

% 将物质点均匀地洒在球面上
clc;clear
r = 1;
X = rand(1000,3);
nt = 500;
c = 100; % 排斥常数
m = 1; % 电荷
dx = sqrt(4*pi*r^2)/(sqrt(size(X,1)));
delta = 3*dx;
bc = 12.0d0 * c / (pi * delta^4);
dt = 10;
massvec = 0.25d0 * dt * dt * (pi*delta^3) * bc / dx;
bforce = zeros(size(X));
velhalfold = zeros(size(X));
velhalf = zeros(size(X));
pforceold = zeros(size(X));
dis = zeros(size(X));
qq = true(size(X,1),1);
for tt = 1:1:nt
    fprintf('tt=%d\n',tt);
    Ltemp = vecnorm(X+dis,2,2);
    if(mod(tt,10)==0)
        figure(4)
        clf
        DT = delaunayTriangulation(X);
        [ccc,v] = convexHull(DT);
        patch('Faces',ccc,'Vertices',X,'FaceColor','cyan')
        axis equal
        axis off
        title(['tt=',num2str(tt)],'fontsize',16)
        pause(0.000001)
    end
    X = (X + dis)./[Ltemp, Ltemp, Ltemp]*r;
    F = zeros(size(X));
    for i = 1:1:size(X,1)
        qq(i) = false;
        C = repmat(X(i,:),size(X,1),1);
        Y = X(qq,:) - C(qq,:);
        L = vecnorm(Y,2,2);
        MAG = L.^3/(c*m*m);
        M = -Y./[MAG,MAG,MAG];
        F(i,:) =  sum(M,1) ;
        qq(i) = true;
    end
    % 粒子沿着球切向受力
    N = sum(F.*X/r,2);
    F = F - [N,N,N].*X/r;
    pforce = F;
    % 动态松弛法
    q = abs(velhalfold)>1e-12;
    cn1 = - sum(sum(dis(q) .* dis(q) .* (pforce(q) / massvec - pforceold(q) / massvec) ./ (dt * velhalfold(q))));
    cn2 = sum(sum(dis.*dis));
    if(abs(cn2)>1e-12)
        if ((cn1 / cn2) > 0.0)
            cn = 2.0 * sqrt(cn1 / cn2);
        else
            cn = 0.0;
        end
    else
        cn = 0.0;
    end
    if (cn > 2.0d0)
        cn = 1.9d0;
    end
    %cn = 1.5;
    if (tt==1)
        velhalf = 1.0d0 * dt / massvec * (pforce + bforce) / 2.0;
    else
        velhalf = ((2.0d0 - cn * dt) * velhalfold + 2.0 * dt / massvec * (pforce + bforce)) / (2.0 + cn * dt);
    end
    % 时间积分
    vel = 0.5d0 * (velhalfold + velhalf);
    dis = dis + velhalf * dt;
    velhalfold = velhalf;
    pforceold = pforce;
end
% 最后X就是你需要的
figure(5)
clf
scatter3(X(:,1),X(:,2),X(:,3),10,'r','filled')
axis equal
axis off

效果：

如有帮助还望题主给个采纳支持一下答主答题哟，谢谢啦