目前有找到任何简单方便的公式得到球面上均匀分布的n点的坐标，在实际工作中可以使用如下方法生成近似均匀分布的点坐标（matlab代码）

%% **************************** 在单位球面上生成均匀分布的n个点 ***********************%
% 参考文献：
% 1. https://www.cnblogs.com/cofludy/p/5894270.html
% 2. https://zhidao.baidu.com/question/2266363508196979468.html
clc;
close all;
clear all;
fclose('all');
format shortG;

nums_points=100; % 点数量
nums_interas=200; % 迭代次数
monitor = 1; %迭代过程监控开关，0：off；1：on

r = get_points_from_sphere(nums_points, nums_interas, monitor);

figure
plot3(r(:,1),r(:,2),r(:,3),'.');hold on;%画结果
[xx,yy,zz]=sphere(50); 
h2=surf(xx,yy,zz); %画一个单位球做参考
set(h2,'edgecolor','none','facecolor','r','facealpha',0.7);
axis equal;
axis([-1 1 -1 1 -1 1]);
hold off;

function r = get_points_from_sphere(nums_points, nums_interas, monitor)
% 在单位球面上生成均匀分布的nums_points个点
% 输入：
% nums_points  点数量
% nums_interas 迭代次数
% monitor  迭代过程监控开关，0：off；1：on
% 输出：
% nums_points个点的坐标

% 参考文献：
% 1. https://www.cnblogs.com/cofludy/p/5894270.html
% 2. https://zhidao.baidu.com/question/2266363508196979468.html
a=rand(nums_points, 1)*2*pi;  %根据随机求面均匀分布，先生成一个初始状态
b=asin(rand(nums_points, 1)*2-1);
r=[cos(a).*cos(b), sin(a).*cos(b), sin(b)];
v=zeros(size(r));
G=1e-2;%斥力常数，试验这个值比较不错
Distance_limit = sqrt(4 * pi / nums_points / pi) * 0.01;
if monitor ==1
    figure
    hold on;
    title('收敛过程'); % Set title
    xlabel("迭代次数", 'fontsize', 10); % Set X axis title
    ylabel("各点间距之和", 'fontsize', 10); % Set Y axis title
end
for ii=1 : nums_interas % 模拟200步，一般已经收敛，其实可以在之前退出    
    dd=zeros(3,nums_points,nums_points); %各点间的矢量差
    
    for m=1:nums_points-1
          for n=m+1:nums_points
              dd(:,m,n)=(r(m,:)-r(n,:))';
              dd(:,n,m)=-dd(:,m,n);
          end
    end 
    L=sqrt(sum(dd.^2,1)); % 各点间的距离
    
    if monitor ==1
        plot(ii, sum(sum(L))/2,'+','color', 'r');
        pause(0.001);
        drawnow;
    end
    L(L<Distance_limit)=Distance_limit; % 距离过小限定
    F=sum(dd./repmat(L.^3,[3 1 1]),3)'; % 计算合力
    Fr=r.*repmat(dot(F,r,2),[1 3]); % 计算合力径向分量
    Fv=F-Fr; % 切向分量
    r=r+v;  % 更新坐标
    r=r./repmat(sqrt(sum(r.^2,2)),[1 3]);
    v=v+G*Fv; % 更新速度
end

end