CODEWITHSUNDEEP
pythonintersectioncurveplane
Mike Cy
Mike Cy

Reputation: 41

Helix-plane intersection

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.

Upvotes: 0

Views: 586

Answers (1)

bombadilhom
bombadilhom

Reputation: 131

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.

Upvotes: 1

Related Questions