首先，您需要在圆的平面中定义两个基本向量。
第一个是正交于法线n = p2-p1的任意矢量
选择具有最大幅值的法线分量和具有第二幅值的分量。
交换它们的值，取最大值的反，并使第三个分量为零（注意，结果与法线的点积为零，因此它们是正交的）
例如，如果n.y最大，n.z次之，则使

v = (0, n.z, -n.y)

然后使用矢量积计算第二个基矢量

u = n x v

规范化向量v和u。使用向量形式上的中心点p1圈出点：

f(rho) = p1 + r * v * cos(rho) + r * u * sin(rho)

或在组件中：

f.x = p1.x + r * v.x * cos(rho) + r * u.x * sin(rho)
and so on

假设我们有一个向量n，我们想找到一个以圆心p1为圆心、半径为r、与n正交的点组成的圆。

p1 = np.array([-21.03181359,   4.54876345,  19.26943601])
n = np.array([-0.06592715,  0.00713031, -0.26809672])
n = n / np.linalg.norm(n) # normalise n
r = 0.5

x = np.array([1,0,0]).astype(np.float64) # take a random vector of magnitude 1
x -= x.dot(n) * n / np.linalg.norm(n)**2  # make it orthogonal to n
x /= np.linalg.norm(x)  # normalize

# find first point on circle (x1). 
# currently it has magnitude of 1, so we multiply it by the r
x1 = p1 + (x*r)

# vector from lumen centre to first circle point
p1x1 = x1 - p1

def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

# rotate the vector p1x1 around the axis n with angle theta
circle = []
for theta in range(0,360,6):
    circle_i = np.dot(rotation_matrix(n, np.deg2rad(theta)), p1x1)
    circle.append(circle_i+p1)

ax = axes3d.Axes3D(plt.figure(figsize=(10,10)))
ax.scatter3D(*np.array(circle).T, s=10, c='red')
ax.scatter3D(*p1.T, s=10, c='black')
ax.set_xlabel('X', size=40)
ax.set_ylabel('Y', size=40)
ax.set_zlabel('Z', size=40)

ax.set_xlim(-19,-22)
ax.set_ylim(2,5)
ax.set_zlim(18,21)