螺旋平面相交答案

【问题标题】：Helix-plane intersection螺旋平面相交
【发布时间】：2018-12-28 17:18:14
【问题描述】：

如何使用（修改）Python 代码找到螺旋线（x=Rcos(t)、y=Rsin(t)、z=a*t）与平面（n - 平面的法线向量和 p0）的交点- 平面上的点）？谢谢。在帖子'3D Line-Plane Intersection'中有答案如何为由两点定义的线做这样的事情，但我需要螺旋线的解决方案。

【问题讨论】：

注意解可能不是唯一的。

标签： python intersection curve plane

【解决方案1】：

您需要解方程 (h(t)-p0).n = 0，其中 h(t) 是您的螺旋。

这个方程没有简单的解析解，但你可以用数值求解，例如 scipy：

import numpy as np
from scipy import optimize

n = np.array([nx, ny, nz])
p0 = np.array([p0x, p0y, p0z])

def h(t):
    return np.array([R*np.cos(t), R*np.sin(t), a*t])

res = optimize.minimize_scalar(lambda t: np.dot(h(t) - p0, n))

print(res.x)

如果你没有scipy/numpy，在这种特定情况下实现Newton method相对容易（我们可以解析计算h(t)的导数）。纯python版本：

from math import cos, sin

n = [nx, ny, nz]
p0 = [p0x, p0y, p0z]

def dot(a, b):
    return sum([x*y for x, y in zip(a, b)])

def h(t):
    return [R*cos(t), R*sin(t), a*t]

def hp(t): # the derivative of h
    return [-R*sin(t), R*cos(t), a]

def find_root_newton(x, f, fp, epsilon=1e-5):
    xn = x + 2*epsilon
    while(abs(xn - x) > epsilon):
        x = xn
        xn = x - f(x)/fp(x)
    return xn

t = find_root_newton(0., lambda t: dot(h(t), n) - dot(p0, n),
                     lambda t: dot(hp(t), n))
print(h(t))

如果螺旋的轴在平面上，它可能会失败（在这种情况下，无论如何你的问题定义不好），而且效率不高。

【讨论】：

感谢您的帮助。会尝试。我没有提到它应该与没有 SciPy 的 Abaqus 中的 Python 一起使用。我将首先安装 SciPy，然后看看它是否可以在 Abaqus 中运行。还有其他方法吗？
您好，bombadilhom，非常感谢您的帮助。您说“Numpy 也不是必需的”，但在代码中可以看到“np”。那么，这意味着什么？此外，当我运行代码时，我得到“NameError：未定义全局名称'h'”。这是什么？
我的意思是，虽然我编写的代码使用了 numpy，但它很容易不用。 'h' 在代码的第一部分中定义。我编辑了帖子以编写一个不使用 numpy 的自包含版本。
代码完全按照它应该的方式工作。非常感谢您的帮助。
经过额外验证，发现计算不正确：附图中，计算点在平面上（Y-X视图），而不是在螺旋上（Y-Z视图）。如何更改代码以获得正确的结果？谢谢你。（不知道如何附加 png 文件）。