Keeping the arc length constant between points in a spiral

Question

I'm making a visualization of points in a logarithmic spiral and want to keep the arc length between points (image particles) constant. I read that in an Archemedian spiral arc length is s = r * theta. I've been fiddling with that and can't quite understand how to implement it (replace theta with s/r where s is an arbitrary constant and increment only r?), and then how to apply it to the case of a logarithmic spiral.

Here are my two test case spirals below, with images. I'd appreciate a nudge in the right direction on how to get these image particles evenly spaced throughout each spiral form. Kind thanks!

An Archemedian spiral, with image:

// ARCHEMEDIAN SPIRAL
var r = 1;
var theta = 1;
for (var p = 0; p < particleCount; p++) {
    var pX = r * Math.cos(theta);
    var pY = r * Math.sin(theta);
    theta += 1;
    r += 10;
    // make an image particle at pX, pY

A logarithmic spiral, with image:

var r = 1;
var theta = 1;
var a = 100; // a constant
var b = 0.2; // another constant
for (var p = 0; p < particleCount; p++) {
    var pX = a * Math.cos(theta) * Math.exp(b * theta);
    var pY = a * Math.sin(theta) * Math.exp(b * theta);
    theta += 0.1;
    // make an image particle at pX, pY

Answer 1

Logarithmic spiral

Approximate solution

The logarithmic spiral is self-similar. As a consequence of that, the tangents make the same angle with the position vector in each point of the spiral. You can use that to draw a reasonable approximation of the spiral simply by taking unit steps in the direction of the current tangent.

Perhaps it's easiest to look at this in terms of complex numbers:

$$ z(\theta)=ae^{(i+b)\theta}=a\left(e^{i+b}\right)^\theta\\ \frac{\mathrm dz}{\mathrm d\theta} = ae^{(i+b)\theta}(i+b) $$

So the direction vector is the position vector rotated by $(i+b)$, which means the angle between them is

$$\alpha=\arctan\frac{1}{b}=\operatorname{arccot}b $$

So the above approach would be something like this (written in sage code):

def logspiral1(s):
    hypot = sqrt(1 + b^2)
    sinA = 1/hypot
    cosA = b/hypot
    x = 1
    y = 0
    for c in range(100):
        yield (x, y)
        f = s/sqrt(x^2 + y^2)
        dx = f*(cosA*x - sinA*y)
        dy = f*(sinA*x + cosA*y)
        x += dx
        y += dy

Looks like a log spiral if taken all by it self, but since errors accumulate it will diverge considerably from the parametric equation.

Exact solution

To avoid that, you should probably compute the new $\theta$ in terms of the direction you want to take. Doing so accurately would mean integrating curve length.

$$\int_{\theta_1}^{\theta_2}\left\lvert ae^{(i+b)\theta}(i+b)\right\rvert\,\mathrm d\theta$$

According to Wolfram Alpha this can be computed using

$$\int\left\lvert ae^{(i+b)\theta}(i+b)\right\rvert\;\mathrm d\theta= \frac{\lvert a(b+i)\rvert}{b}e^{b\theta}$$

so you want

$$ s=\frac{\lvert a(b+i)\rvert}{b}\left(e^{b\theta_2}-e^{b\theta_1}\right) \\ e^{b\theta_2}=\frac{sb}{\lvert a(b+i)\rvert}+e^{b\theta_1} \\ \theta_2=\frac1b\log\left(\frac{sb}{\lvert a(b+i)\rvert}+e^{b\theta_1}\right) $$

def logspiral2(s):
    theta = 0
    h = s*b/(a*sqrt(1 + b^2))
    for c in range(100):
        yield (a*cos(theta)*exp(b*theta), a*sin(theta)*exp(b*theta))
        theta = log(h + exp(b*theta))/b

The first, approximate approach is drawn in red in the following figure, the second exact one in green. The continuous line behind the green points is a simple parametric plot of the exact solution.

Archimedean spiral

Attempted exact solution

So let's do the same for the Archimedean spiral, again using Wolfram Alpha.

$$z(\theta)=r\theta e^{i\theta} \\ \frac{\mathrm dz}{\mathrm d\theta} = r(1+i\theta)e^{i\theta} \\ \int\left\lvert r(1+i\theta)e^{i\theta}\right\rvert\;\mathrm d\theta = \int r\sqrt{1+\theta^2}\;\mathrm d\theta =\frac r2\left(\theta\sqrt{1+\theta^2}+\operatorname{arsinh}\theta\right)$$

That last term is a real beast: it includes $\theta$ both inside and outside the transcendental function $\operatorname{arsinh}$ (the inverse sinus hyperbolicus), which means we probably won't be able to solve this expression for $\theta$ except perhaps numerically.

Approximate solution

You could again do this thing we did above, take unit steps in tangent direction. Since the tangent direction angle depends on $\theta$, you'd have to compute that:

def archspiral1(s):
    x = 0
    y = 0
    for c in range(100):
        yield (x, y)
        theta = float(sqrt(x^2 + y^2)/r)
        f = float(s/sqrt(1 + theta^2))
        dx = f*(1*x - theta*y)
        dy = f*(theta*x + 1*y)
        if theta == 0: dx = s
        x += dx
        y += dy

Exact curve but approximate distances

Or you could consider the change in $\theta$ induced by such a step in tangent direction. If you look at the derivative, then you notice that the radial component of that distance is $r$ while the tangential (to the circle) component is $r\theta$. So you have something like

$$s = x\sqrt{1+\theta^2}\qquad r\left(\theta_2-\theta_1\right) = x \\ \theta_2=\theta_1+\frac s{r\sqrt{1+\theta_1^2}}$$

def archspiral2(s):
    theta = s/2
    for c in range(100):
        yield (r*theta*cos(theta), r*theta*sin(theta))
        theta += s/(r*sqrt(1 + theta^2))