Helix-plane intersection

Question

How can I use (modify) Python code to find intersection of helix (x=Rcos(t), y=Rsin(t), z=a*t) with plane (n - normal vector of the plane and p0 - point on plane)? Thanks. In post '3D Line-Plane Intersection' there are answers how to do such thing for line defined by two points but I need solution for helix.

Note that the solution may not be unique. – Richard Dec 28 '18 at 18:09 — Richard, Dec 28 '18 at 18:09

bombadilhom · Answer 1 · 2018-12-30T12:53:45.123

1

You'll need to solve the equation (h(t)-p0).n = 0, where h(t) is your helix.

The equation does not admit a easy analytic solution, but you can solve it numerically, with scipy for example:

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)

If you don't have scipy/numpy, it's relatively easy to implement the Newton method in this specific situation (we can analytically compute the derivative of h(t)). Pure python version:

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))

It can fail if the axis of the helix is in the plane (in which case your problem is badly defined anyhow), and it's not really efficient.

edited Dec 30 '18 at 12:53

answered Dec 28 '18 at 17:32

bombadilhom

131
1
4

Thank you for the help. Will try. I didn't mention that it should works with Python in Abaqus which has no SciPy. I'll first install SciPy and will see if it will work in Abaqus. Is there another way around? – Mike Cy Dec 29 '18 at 01:56
Hi bombadilhom, Thanks a lot for your help. You said "Numpy is not necessary either" but in code one can see 'np'. So, what does it mean? Also, when I run code, I've got "NameError: global name 'h' is not defined". What is it? – Mike Cy Dec 30 '18 at 02:48
I meant that while the code I wrote used numpy it was easy to do without. 'h' was defined in the first part of the code. I edited the post to write a self contained version that does not use numpy. – bombadilhom Dec 30 '18 at 12:56
The code works exactly as it should. Thank you very much for the help. – Mike Cy Dec 30 '18 at 16:52
After additional verification, I found that the calculations is incorrect: in the attached picture, the computed point lies on the plane (Y-X View), but not on the helix (Y-Z View). How can I change the code to get the correct result? Thank you. (Don't know how to attach png-file). – Mike Cy Jan 02 '19 at 22:45
After all everything is OK: in test example I had negative sign for Z and instead z=-a*t , I've used z=a*t. – Mike Cy Jan 03 '19 at 02:09

Helix-plane intersection

1 Answers1