在平面上均匀分布的点构成正三角形网格，在平面上划出正二十面体的表面。若将正二十面体叠好，则原均匀分布于平面上的点会均匀分布在正二十面体表面。将正二十面体的表面保面积地映到球面上，则得到球面上近似均匀分布的点。
详见 https://doi.org/10.48550/arXiv.1607.04590
使用python实现的代码如下：
def EqualAreaIcosahedral(m, n):
"""
D.P. Hardin, T. Michaels and E.B. Saff
A Comparison of Popular Point Configurations on S2.
October 25, 2016
Parameters:
m - integer
n - integer
Returns:
a list of 3-element list, every 3-element list means the 3D cartesian
coordinate of a point
"""
import math
def rotate(r, n, phi):
# Rodrigues Formula
# r is a 3D vector, n is a unit vector describing an axis of rotation
# about which r rotates by an angle phi according to the LEFT hand rule
# return the rotated vector
cosphi = math.cos(phi)
sinphi = math.sin(phi)
nr = n[0]*r[0] + n[1]*r[1] + n[2]*r[2]
x = r[0]*cosphi + n[0]*nr*(1-cosphi) + (r[1]*n[2]-r[2]*n[1])*sinphi
y = r[1]*cosphi + n[1]*nr*(1-cosphi) + (r[2]*n[0]-r[0]*n[2])*sinphi
z = r[2]*cosphi + n[2]*nr*(1-cosphi) + (r[0]*n[1]-r[1]*n[0])*sinphi
# print(x*x+y*y+z*z)
return (x,y,z)
s_3 = math.sqrt(3)
M = max(m, n)
N = min(m, n)
assert N >= 0 and M > 0
# the length of icosahedral's edge
L1 = 2*math.sqrt(math.pi)/math.sqrt(math.sqrt(75))
L2 = M*M+N*N+M*N
# the greatest common divisor of M & N
for i in range(N):
if M%(i+1) == 0:
gcd = i+1
if N == 0:
gcd = M
# whether there is a point at the center of icosahedral's face
flag_center = (M-N)%3 == 0
# points on the first 2/3 angular bisectors
L60 = [[], [], []]
# cartesian coordinates of points
p1 = []
for i in range(-N, M+1):
for j in range(M+N+1):
# i & j - skew coordinates of points
# i & j in diamond contain the equilateral triangle
if flag_center and 3*i == (M-N) and 3*j == (M+2*N):
continue
x = i + j/2
y = j*s_3/2
# 3 angular bisector
if j>0 and i*(M+2*N) == (M-N)*j and M+2*N > 3*j:
L60[0].append((x,y))
continue
if i<M and (2*M+N)*(j-N) == (M-N)*(M-i) and 2*M+N > 3*(M-i):
L60[1].append((x,y))
continue
if j<M+N and (M+2*N)*(M+N-j) == (i+N)*(2*M+N) and (M+2*N) > 3*(i+N):
L60[2].append((x,y))
continue
# check whether points in the equilateral triangle
x1 = (M+N)*x*s_3 - (M-N)*y
x2 = N*x*s_3 - (2*M+N)*y
x3 = M*x*s_3 + (M+2*N)*y
if x1 > 0 and x2 < 0 and x3 < L2*s_3:
p1.append((x,y))
# cartesian coordinates of points, while rotate the equilateral triangle so that
# its edge coincide with the x-axis and with length 1
p2 = []
for p in p1:
x = p[0]*(M+N/2) + p[1]*N*s_3/2
y = p[1]*(M+N/2) - p[0]*s_3*N/2
p2.append((x/L2, y/L2))
# h & w in reference paper
hw = [[], [], []]
for p in p2:
if p[1] < p[0]/s_3 and p[0] + p[1]*s_3 <= 1:
h = s_3/6 - p[1]
w = p[0] - 0.5
hw[0].append((h*L1,w*L1))
elif p[0] < 0.5:
h = (p[1] - p[0]*s_3 + 1/s_3)/2
w = (1 - p[0] - p[1]*s_3)/2
hw[1].append((h*L1,w*L1))
else:
h = (p[0]*s_3 + p[1])/2 - 1/s_3
w = (p[1]*s_3 - p[0])/2
hw[2].append((h*L1,w*L1))
# lambda, psi in reference paper
lp = [[], [], []]
for i in range(3):
for p in hw[i]:
ps = p[0]*p[0]*s_3/2 + math.pi/6
la = math.atan(math.sin(p[0]*p[1]/2)/(math.cos(p[0]*p[1]/2)-2*math.cos(ps)/s_3))
lp[i].append((la,ps))
for p in L60[0]:
y = (p[1]*(M+N/2) - p[0]*s_3*N/2) / L2
h = (s_3/6 - y) * L1
ps = h*h*s_3/2 + math.pi/6
lp[i].append((math.pi/3,ps))
# theta, phi - spherical coordinates of points
tp = []
for i in range(3):
for p in lp[0]:
cosPhi = math.cos(p[1])*math.sin(p[0])*2/s_3
# print(math.cos(p[0]), math.pi/p[0], math.pi/p[1], math.sqrt(1-cosPhi**2),
# math.cos(p[0])*math.cos(p[1])*2/math.sqrt(3*(1-cosPhi**2)))
th = math.acos(math.cos(p[0])*math.cos(p[1])*2/math.sqrt(3*(1-cosPhi**2)))
tp.append((th, p[0]+i*math.pi*2/3))
# 3D cartesian coordinates of vertex and points on the 1st ~ 20th face
xyz = [[] for i in range(21)]
# points on the first face
for p in tp:
x = math.sin(p[0]) * math.cos(p[1])
y = math.sin(p[0]) * math.sin(p[1])
z = math.cos(p[0])
xyz[1].append((x,y,z))
# north pole
if flag_center:
xyz[1].append((0,0,1))
# rotate and reflex to obtain the else ponits
cosc = 1/(s_3*math.tan(math.pi/5))
sinc = math.sqrt(1-cosc**2)
for p in xyz[1]:
for i in range(4):
xyz[i+2].append(rotate(p,(-sinc,0,cosc),math.pi*(i+1)*0.4))
xyz[6].append(rotate(p,(sinc/2,s_3*sinc/2,cosc),-math.pi*0.4))
xyz[11].append(rotate(p,(sinc/2,s_3*sinc/2,cosc),-math.pi*0.8))
xyz[16].append((-p[0],-p[1],-p[2]))
for p6, p11, p18 in zip(xyz[6], xyz[11], xyz[16]):
for i in range(4):
xyz[7+i].append(rotate(p6,(-sinc,0,cosc),math.pi*(i+1)*0.4))
xyz[12+i].append(rotate(p11,(-sinc,0,cosc),math.pi*(i+1)*0.4))
xyz[17+i].append(rotate(p18,(-sinc,0,cosc),math.pi*(i+1)*0.4))
# vertex
xyz[0] += [(-sinc,0,cosc), (sinc/2,s_3*sinc/2,cosc)]
for a in [math.pi*(i+1)*0.4 for i in range(4)]:
xyz[0].append(rotate((sinc/2,s_3*sinc/2,cosc),(-sinc,0,cosc),a))
for p in xyz[0][:6]:
xyz[0].append((-p[0],-p[1],-p[2]))
# all ponits
xyz.append([])
for ps in xyz[:-1]:
xyz[-1] += ps
# if there are points on edges
if gcd > 1:
xyze = []
hwe = [(s_3*L1/6,(-0.5+(i+1)/gcd)*L1) for i in range(gcd-1)]
lpe = []
for p in hwe:
ps = p[0]*p[0]*s_3/2 + math.pi/6
la = math.atan(math.sin(p[0]*p[1]/2)/(math.cos(p[0]*p[1]/2)-2*math.cos(ps)/s_3))
lpe.append((la,ps))
tpe = []
for p in lpe:
cosPhi = math.cos(p[1])*math.sin(p[0])*2/s_3
th = math.acos(math.cos(p[0])*math.cos(p[1])*2/math.sqrt(3*(1-cosPhi**2)))
tpe.append((th, p[0]))
tpe.append((th, p[0]+math.pi*4/3))
for p in tpe:
x = math.sin(p[0]) * math.cos(p[1])
y = math.sin(p[0]) * math.sin(p[1])
z = math.cos(p[0])
xyze.append((x,y,z))
xyze.append((-x,-y,-z))
xyze.append(rotate((x,y,z),(sinc/2,s_3*sinc/2,cosc),-math.pi*0.4))
for p in xyze[:3*len(tpe)]:
for a in [math.pi*(i+1)*0.4 for i in range(4)]:
xyze.append(rotate(p,(-sinc,0,cosc),a))
xyz[-1] += xyze
return xyz[-1]