9
$\begingroup$

I am looking for a way to place points equidistantly along an Archimedes spiral according to arch-length (or an approximation) given the following parameters:

Max Radius, Fixed distance between the points, Number of points

I have been lurking around on this site for a few days and have found a lot of great advice but am struggling with the syntax of some of the proposed responses (not a native coder but I have had small exposure to Python and Matlab).

This example 1 seems to be exactly what I am looking for but I am just struggling with the code, it is not clear to me what variables are used or how the program executes.

Example 2 and example 3 were also helpful but I am definitely missing something when it comes to solving the equation numerically as the resulting spiral does not have equal spacing.

My goal is to use a spreadsheet (MS Excel) to drive a solid modeling program to generate a hole pattern per the parameters above.

Cheers!

$\endgroup$

3 Answers 3

8
$\begingroup$

I've taken a different approach, which, while requiring an extra step is ultimately a better approximation to an equi-spaced curve. In the complex plane, we know that the arc length is given by

$$s=\int_0^{\theta}|\dot z|d\theta$$

where $\dot z$ means differentiation w.r.t. $\theta$ in this case. The for range of $[\theta:\theta+\Delta \theta]$ we can say

$$\Delta s=\int_{\theta}^{\theta+\Delta \theta}|\dot z|d\theta\approx \frac{|\dot z|_n+|\dot z|_{n+1}}{2}\Delta \theta$$

So, the extra step I mentioned consists of building successive $\theta_n$ from $\theta_{n-1}$. For Archimedes spiral $|\dot z|=\sqrt{1+\theta^2}$ and we can show that

$$\Delta \theta_n\approx\frac{\Delta s}{\sqrt{1+\theta_n^2}}$$

The figure below shows a comparison of the uniform $\theta$ (left) and uniform $s$ (right). I kept the derivation very general up to the last minute because it applies to any curve, obviously. Archimedes Spirals.

$\endgroup$
2
  • $\begingroup$ This answer will be extremely useful in my work developing algorithms for path planning of a polar robot. $\endgroup$
    – Ravenex
    Aug 4, 2017 at 15:06
  • $\begingroup$ @Ravenex Thanks, you'll find a more recent, and more general discussion here: math.stackexchange.com/questions/2375940/…. $\endgroup$ Aug 4, 2017 at 15:26
3
$\begingroup$

If the polar equation of an Archimedean spiral is given by: $$ \rho = \theta $$ then its parametric equation is $(\theta\cos\theta,\theta\sin\theta)$ and the arc length between $0$ and $\theta_f$ is given by: $$ L= \int_{0}^{\theta_f}\sqrt{1+\theta^2}\,d\theta = \frac{1}{2}\left(\theta_f \sqrt{1+\theta_f^2}+\text{arcsinh}(\theta_f)\right)\approx \theta_f \sqrt{1+\frac{\theta_f^2}{4}}$$ so a good way to place almost-equispaced point on such Archimedean spiral is to take the $N$-th point at: $$ \theta_N = \sqrt{2}\sqrt{-1+\sqrt{1+k^2 N^2}}. $$

$\endgroup$
2
  • 1
    $\begingroup$ As a simpler alternative, one may consider $\theta_N=\sqrt{N/N_{\rm points}}\theta_{\rm MAX}$, which is fairly accurate except near to the origin. $\endgroup$ Jul 24, 2015 at 10:02
  • 4
    $\begingroup$ Uhm, what is k? $\endgroup$
    – Ghanima
    Aug 24, 2017 at 16:17
-1
$\begingroup$

A python implementation to solve for evenly spaced points around an Archimedean spiral. Most of the solutions I came across assume $a=0$ so here is my solution which does not ignore $a$.

Use a solver to iteratively find $\theta_b$ that satisfies:

$\text{arc length} = s = \int_{\theta_a}^{\theta_b} \sqrt{(\frac{dr}{d\theta})^2 + r^2}$ when $r = a+b\theta$
https://en.wikipedia.org/wiki/Arc_length#Other_coordinate_systems

$\therefore s = \int_{\theta_a}^{\theta_b} \sqrt{b^2 + (a+b\theta)^2}$

Using Sage to integrate:
$s = \frac{b^{2} \operatorname{arsinh}\left(\frac{b \theta_{b} + a}{b}\right) + \sqrt{b^{2} \theta_{b}^{2} + 2 \, a b \theta_{b} + a^{2} + b^{2}} {\left(b \theta_{b} + a\right)}}{2 \, b} -\frac{b^{2} \operatorname{arsinh}\left(\frac{b \theta_{a} + a}{b}\right) + \sqrt{b^{2} \theta_{a}^{2} + 2 \, a b \theta_{a} + a^{2} + b^{2}} {\left(b \theta_{a} + a\right)}}{2 \, b} $

Scipy's fsolve searches for roots where func(x) = 0. So subtract $s$ (the requested arc length) from both sides and solve for the $\theta_b$ where $\text{func}(\theta_b) \approx 0$.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve


def arc_len(theta_b: float, theta_a: float, a: float, b: float, arclen: float):
    """calculated arc_len between theta_a and theta_b (where theta_b > theta_a) minus requested arclen

    Args:
        theta_b (float): upper bound for arc
        theta_a (float): lower bound for arc
        a (float): from the definition of an archimedean spiral r = a + b * theta
        b (float): from the definition of an archimedean spiral r = a + b * theta
        arclen (float): the goal arc length

    Returns:
        float: arclen from theta_a to theta_b minus requested arc length

    Arc length definition:
    https://en.wikipedia.org/wiki/Arc_length#Other_coordinate_systems

    using sagemath: https://www.sagemath.org/
    sage: arclen, a, b, t, theta_a, theta_b = var('arclen a b t theta_a theta_b')
    sage: integral((b^2 + (a + b*t)^2)^0.5, t, theta_a,theta_b, assume(a>0,b>0,theta_b-theta_a>0)).simplify()

    .. math::
        -\frac{b^{2} \operatorname{arsinh}\left(\frac{b \theta_{a} + a}{b}\right) + \sqrt{b^{2} \theta_{a}^{2} + 2 \, a b \theta_{a} + a^{2} + b^{2}} {\left(b \theta_{a} + a\right)}}{2 \, b} + \frac{b^{2} \operatorname{arsinh}\left(\frac{b \theta_{b} + a}{b}\right) + \sqrt{b^{2} \theta_{b}^{2} + 2 \, a b \theta_{b} + a^{2} + b^{2}} {\left(b \theta_{b} + a\right)}}{2 \, b}

    """

    return (1 / (2 * b)) * (
        (
            b**2 * np.arcsinh((b * theta_b + a) / b)
            + np.sqrt(b**2 * theta_b**2 + 2 * a * b * theta_b + a**2 + b**2)
            * (b * theta_b + a)
        )
        - (
            b**2 * np.arcsinh((b * theta_a + a) / b)
            + np.sqrt(b**2 * theta_a**2 + 2 * a * b * theta_a + a**2 + b**2)
            * (b * theta_a + a)
        )
    ) - arclen


def main():
    # archimedean spiral definition r = a + b*theta
    a = 3
    b = 1 / (2 * np.pi)
    turns = 3

    # requested arc length
    arclen = 7

    theta = 0
    max_theta = turns * 2 * np.pi

    thetas = []

    while theta < max_theta:
        theta = fsolve(arc_len, [theta], args=(theta, a, b, arclen,))[0]
        thetas.append(theta)

    # ignore the last calculated point (Outside of the requested range)
    thetas = np.array(thetas[:-1])

    ################ PLOTTING ################
    ### polar ###
    fig = plt.figure(figsize=(6,3))
    ax0 = plt.subplot(121, polar=True, aspect="equal")

    # calculated points
    r_div = a + b * thetas
    ax0.scatter(thetas, r_div, color="#a65628", zorder=2)

    # spiral line
    theta = np.arange(0, turns * 2 * np.pi + 0.001, np.pi / 32)
    r = a + b * theta
    ax0.plot(theta, r, color="#377eb8", zorder=1)

    ### cartesian ###
    ax1 = plt.subplot(122, aspect="equal")

    x = r_div * np.cos(thetas)
    y = r_div * np.sin(thetas)

    # calculated points
    ax1.scatter(x, y, color="#a65628", zorder=2)

    x = r * np.cos(theta)
    y = r * np.sin(theta)

    # spiral line
    ax1.plot(x, y, color="#377eb8", zorder=1)

    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Plotted in both polar and Cartesian, the brown points here are spaced evenly (with an arc length = 7) around the spiral $r = 3 + \frac{1}{2\pi}\theta$

figure created by the above code

$\endgroup$
1
  • $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Apr 6 at 8:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.