How do I calculate this loop spline given the length, angle and horizontal offset?

Question

I'm developing a formula to calculate a loop spline from a length, angle and horizontal offset. I can successfully calculate the loop from the first two parameters, but taking the horizontal offset into account is over my head. I can simply add the horizontal offset, but that will make the spline longer and no longer match with the required "length" parameter.

Assuming $\theta = angle \neq 0$ and $x \in [0,1]$, the formula I use is this

$$ \begin{align} r &= \dfrac{\dfrac{length}{2\pi}}{\dfrac{\theta}{2\pi}} = \dfrac{length}{\theta} & (1)\\[2ex] p_x &= cos(x * \theta)*r \\ p_y &= sin(x * \theta)*r & (2)\\ p_z &= x * h_{offset} \end{align} $$

(1) Radius of a circle $r = \dfrac{C}{2\pi}$
(2) Calculating the $xyz$ coordinates of the circle in 3D space, gradually offsetting the z-coordinate with the progress of the loop spline to finally meet the horizontal offset

(Disclaimer: My use of mathematical notation is probably far off from how its used properly.)

---

This gives a nice loop spline, but unfortunately, given $h_{offset} \ne 0$, the resulting spline does not match the $length$ requirement. Example below:

---

The spline also needs to be smoothed in the x-component (as currently there's no soft transition when traveling in z-direction and entering the loop, that is coming from the bottom of the grid from the pictures perspective).

This is a function I use to smooth the $h_{offset}$ in each step, resulting in a smooth transition.

$$ x \in [0,1], \ f(x) = \begin{cases} \dfrac{4x^2}{2} & \text{if } x < \dfrac{1}{2} \\ 1-\dfrac{4(1-x)^2}{2} & \text{otherwise} \end{cases} \ \ \ \ (3) $$

$$ p_z = f(x) * h_{offset} \ \ \ (4) $$

(3) I hope I translated the formula correctly, here's the Python version of it (4) New formular for computing the Z component

I think that applying $f()$ to $p_x$ should also have an effect on the length of the spline, but I'm not sure how.

---

Ultimately, my question is: How can I calculate the radius so that in the end the length of the spline matches the $length$ parameter for any given $hoffset$?

I'd be happy to correctly calculate the radius for the version without smoothing out the loop with $f(x)$ as a starter.

---

Edit#1: I've taken samples of the actual length of the spline. From the plot, it appears the length that the horizontal offset adds is marginal at first and then smoothly transitions into a linear function...?

---

Update: June 18th

So basically all I need is a way to calculate the formula that gives the above graph in order to subtract it from the input $lenght$ parameter to counteract the effect of the $h_{offset}$.

I've found this formula to come relatively close to what you can see above in terms of the shape, but I haven't been able to find the exact formula, and using trial & error isn't a very sophisticated approach to the problem.

$$ f(x)=\dfrac{x^2}{x+1} $$

---

By trial & error I figured the below formula matches the above graph very well. But 1) it's not exact and 2) I still can't derive the right formula for arbitrary input parameters.

Answer 1

If we have a parametrized curve $\vec{p}(t) = \left(x(t), y(t), z(t)\right)$, its length (from $t=0$ to $t=1$) is $$L = \int_{0}^{1} \sqrt{ \left.\left(\frac{d x(\tau)}{d \tau}\right)^2\right\rvert_{\tau=t} + \left.\left(\frac{d y(\tau)}{d \tau}\right)^2\right\rvert_{\tau=t} + \left.\left(\frac{d z(\tau)}{d \tau}\right)^2\right\rvert_{\tau=t}} dt$$ This is correct for all three-dimensional parametrized curves.

Here, we have curve $$\vec{p}(t) = \begin{cases} x(t) = r \cos(t \theta) & \, \\ y(t) = r \sin(t \theta) & \,\\ z(t) = 2 h t^2, & t \lt \frac{1}{2} \\ z(t) = h - 2 h (1-t)^2, & t \ge \frac{1}{2} \end{cases}$$ The signs of $\theta$ and $h$ do not affect the curve length for $t=0\dots1$, because their signs only set the orientation (counterclockwise or clockwise) and direction (along the $z$ axis).

If we work out the integral (I use Maple for the heavy lifting), and note that the signs of $\theta$ and $h$ do not affect the curve length, we get (fixed on 2016-07-24): $$L = \begin{cases} \lvert r \theta \rvert, & h = 0 \\ \frac{r^2 \theta^2}{4\lvert h \rvert} \ln\left(\frac{\sqrt{r^2 \theta^2 + 4 h^2} \,+\, 2 \lvert h \rvert}{\lvert r \theta \rvert}\right) \,+\, \frac{1}{2}\sqrt{r^2 \theta^2 + 4 h^2}, & h \ne 0 \end{cases}$$

To find the radius $r$ corresponding to a desired curve length $L_{Desired}$, you'll need to find the solution numerically. I would probably use a binary search, with $\Delta = \frac{L_{Desired}}{\lvert\theta\rvert}$; that is, I'd find the integer $n$, $n \ge 1$, for which $$\begin{cases}L\left(r = n \Delta\right) \le L_{Desired} \\ L\left(r = (n+1) \Delta\right) \gt L_{Desired}\end{cases}$$ and then find the $L(r)=L_{Desired}$ using a simple binary search within $r = [n \Delta \dots (n+1) \Delta]$.

---

The Maple snippet I used to evaluate the (fixed) integral is below.

Note that I assume $h\gt0$ and $\theta\gt0$ here, to simplify the expressions, and add the $\lvert\rvert$'s later.

Also, Niklas R mentioned in a comment that he evaluated only the first half of the integral; this is a smart move (and valid because the curve is symmetric wrt. $t=1/2$), and I also apply it here; the length is calculated as twice the length of the first half.

restart:
interface('showassumed=0'):
assume(t::real, r::real, r>0, h::real, h>0, theta::real, theta>0):

x := t -> r*cos(t*theta);
y := t -> r*sin(t*theta);
z := t -> 2*h*t^2;

L := 2*int( sqrt( diff(x(t),t)^2 + diff(y(t),t)^2 + diff(z(t),t)^2 ), t = 0 .. 1/2);

# with ||, and multiplied by 4*h for simpler form
combine(subs({h=abs(h), theta=abs(theta)}, simplify(L*4*h)));

The first set of lines is just Maple setup and assumptions. The next three are the $x$, $y$, and $z$ functions for $t = 0\dots1/2$. The $L$ line calculates the integral itself. To "prettify" it the way I like it, the final line substitutes back the absolute values of $h$ and $\theta$, and multiplies the expression by $4 h$.

The final output is

r^2 * theta^2 * ln( (2*abs(h) + sqrt(r^2 * theta^2 + 4 * h^2)) / (r * abs(h)) ) + 2 * abs(h) * sqrt(r^2 * theta^2 + 4 * h^2)

If we divide it back by that $4 h$, we get

(r^2 * theta^2) / (4 * abs(h)) * ln( (sqrt(r^2 * theta^2 + 4 * h^2) + 2*abs(h)) / (r * abs(h)) )
+ 1/2 * sqrt(r^2 * theta^2 + 4 * h^2)

which is, as it should be, the (now fixed) value of $L$.

Here is a Python function to calculate L:

import math

def CurveLength(r, theta, h, epsilon = 0.0):

    r = math.fabs(float(r))
    if (r <= epsilon):
        return 0.0

    theta = math.fabs(float(theta))
    if (theta <= epsilon):
        return 0.0

    rtheta = r*theta

    h = math.fabs(float(h))
    if (h <= epsilon):
        return rtheta

    h4 = 4.0*h
    rtheta2 = rtheta*rtheta
    s = math.sqrt(rtheta2 + h*h4)
    return rtheta2*math.log((s + h + h) / rtheta) / h4 + 0.5*s

Here is a Python function that finds the radius corresponding to the desired curve length $L$ (to within caller-specified epsilon) using a binary search:

import math

def FindRadius(L, theta, h, epsilon = 0.0):

    def curve_length(r, theta, h):
        rtheta = r*theta
        rtheta2 = rtheta*rtheta
        h4 = 4.0*h
        s = math.sqrt(rtheta2 + h*h4)
        return rtheta2*math.log((s + 2.0*h) / rtheta) / h4 + 0.5*s

    L = math.fabs(float(L))
    if (L <= epsilon):
        return 0.0

    theta = math.fabs(float(theta))
    if (theta <= 0.0):
        # (epsilon is allowed error in L, so we cannot
        #  compare against it here; we'd need another epsilon
        #  for this.. well, a hard zero will have to do for now.)
        # Cannot achieve L with zero theta.
        # We could instead raise a CallerIsBeingSilly exception.
        return 0.0

    h = math.fabs(float(h))
    if (h <= epsilon):
        return L/theta

    rmin = 0.0
    rmax = L / theta
    rL = curve_length(rmax, theta, h)
    while (rL < L):
        rmin = rmax
        rmax = 2.0 * rmax
        rL = curve_length(rmax, theta, h)

    if (rL - L <= epsilon):
        return rmax

    r = (0.5*rmin) + (0.5*rmax)
    rL = curve_length(r, theta, h)
    while (math.fabs(rL - L) > epsilon):
        if (rL > L):
            rmax = r
        elif (rL < L):
            rmin = r
        else:
            break

        r = (0.5*rmin) + (0.5*rmax)
        rL = curve_length(r, theta, h)

    return r

Note that the latter contains a helper function curve_length, which is the same as the earlier CurveLength, except without the range/sanity checks.

Answer 2

I've found one of my friends being able to provide the answer to this question. It is based on finding the parametrized version of the curve length and solve for the radius.

Let's consider $L := length$ and $h := h_{offset}$

First, we write the three formulas as a function yielding a vector.

$$ p(x) = \begin{bmatrix} cos(\theta x)*r \\ sin(\theta x)*r \\ hx \end{bmatrix} $$

In order to find the length of the curve, we create the derivate $p'(x)$ which represents the velocity vector for the respective point on the curve.

$$ p'(x) = \begin{bmatrix} -sin(\theta x) \theta r \\ cos(\theta x) \theta r \\ h \end{bmatrix} $$

From this function, we derive $p''(x) = |p'(x)|$.

$$ \begin{align} p''(x) &= \sqrt{sin^2(\theta x) \theta^2 r^2 + cos^2(\theta x) \theta^2 r^2 + h^2} \\[0.5em] &= \sqrt{\theta^2 r^2 (sin^2(\theta x) + cos^2(\theta x)) + h^2} \\ &= \sqrt{\theta^2 r^2 + h^2} \end{align} $$

Note how $p''(x)$ is a constant which makes sense since the length of the curve only depends on the $L$ and $h$ parameters. Integrating over $p''(x)$ gives us the length $L$, then solve for $r$

$$ \begin{align} L &= \int_{0}^{1} {p''(x)} dx = p''(x) \\ r &= \sqrt{\dfrac{(L^2-h^2)}{\theta^2}} \end{align} $$

This does not take the smoothing function $f(x)$ into account yet, unfortunately. Following the above pattern with $p_z = f(x) * h$ was too complicated for the both of us. I'd appreciate if anyone has to offer a solution taking $f(x)$ into account as well.

Thanks so far, Max!