Raup’s model

In this tutorial, you will generate coiled shells with Raup’s model (Raup & Michelson, 1965; Raup, 1966).

Raup’s model describes a shell as a trajectory of a generating curve, which mimics the aperture shape, that expands, rotates, and translates along a fixed coiling axis. Three parameters control the geometry:

• whorl expansion rate w_r: how fast the tube widens per revolution

• translation rate t_r: how fast the coil moves down along the coiling axis

• relative distance of the generating curve from the coiling axis d_r

Setup

# Uncomment if needed
# %pip install "ktch[plot]"

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from ktch.coiling import RaupModel, l_r, raup, theta_r

Generate your first shell

The raup function takes the three parameters and returns a sampled surface. We also pass a sampling grid (theta_range, phi_range); omit them to use the defaults.

theta = np.linspace(0.0, 2.0 * np.pi * 4.0, 240)  # four whorls
phi = np.linspace(0.0, 2.0 * np.pi, 60)

X = raup(w_r=1.3, t_r=1.5, d_r=0.2, theta_range=theta, phi_range=phi)
X.shape

(240, 60, 3)

You should see a (n_theta, n_phi, 3) array: n_theta samples along the coil, n_phi samples around the tube, and the last axis holds the (x, y, z) coordinates of each surface point.

Visualize in 3D

Let’s define a small helper that turns a surface array into an interactive plotly figure.

Show definition: plot_surface3d()

Hide definition: plot_surface3d()

def plot_surface3d(X, *, title="", colorscale="Viridis", apex_up=True):
    """Plot a ``(n, n_phi, 3)`` surface array with plotly.

    With ``apex_up=True`` the form is rotated 180° about a horizontal axis so
    the apex points up (the conventional orientation). This is a rotation, not a
    ``z`` reflection, so the coiling direction is preserved.
    """
    X = np.asarray(X, dtype=float)
    if apex_up:
        X = X * np.array([-1.0, 1.0, -1.0])
    x, y, z = X[..., 0], X[..., 1], X[..., 2]
    fig = go.Figure(
        data=[go.Surface(x=x, y=y, z=z, colorscale=colorscale, showscale=False)]
    )
    fig.update_layout(
        title=title,
        width=700,
        height=500,
        scene=dict(aspectmode="data"),
    )
    return fig

fig = plot_surface3d(X, title="Raup's model (w_r=1.3, t_r=1.5, d_r=0.2)")
fig.show()

That’s a coiled shell built from three parameters.

Note

Raup’s model coils along the coiling axis (+z), so the apex sits at the bottom of the raw output. plot_surface3d shows shells apex-up by default (apex_up=True): it rotates the form 180° about a horizontal axis (a rotation, not a z-flip, so the coiling direction is preserved). Pass apex_up=False for the raw orientation.

Explore the parameters

To see what each parameter does, we vary one at a time while holding the other two fixed.

Show definition: compare_raup()

Hide definition: compare_raup()

def compare_raup(values, *, vary, title, **fixed):
    """Generate Raup shells varying one parameter and show them side by side."""
    theta = np.linspace(0.0, 2.0 * np.pi * 4.0, 180)  # four whorls
    phi = np.linspace(0.0, 2.0 * np.pi, 60)  # aperture sampling
    fig = make_subplots(
        rows=1,
        cols=len(values),
        specs=[[{"type": "surface"}] * len(values)],
        subplot_titles=[f"{vary}={v}" for v in values],
        horizontal_spacing=0.01,
    )
    for i, v in enumerate(values, start=1):
        Xi = raup(theta_range=theta, phi_range=phi, **{**fixed, vary: v})
        Xi = Xi * np.array([-1.0, 1.0, -1.0])  # apex-up (see plot_surface3d note)
        fig.add_trace(
            go.Surface(
                x=Xi[..., 0],
                y=Xi[..., 1],
                z=Xi[..., 2],
                showscale=False,
                colorscale="Viridis",
            ),
            row=1,
            col=i,
        )
    scene_names = ["scene"] + [f"scene{k}" for k in range(2, len(values) + 1)]
    fig.update_layout(
        title=title,
        width=900,
        height=350,
        margin=dict(l=0, r=0, t=60, b=0),
        **{name: dict(aspectmode="data") for name in scene_names},
    )
    return fig

Whorl expansion rate w_r

compare_raup([1.1, 2.0, 10.0], vary="w_r", title="Whorl expansion rate",
             t_r=1.0, d_r=0.0).show()

Translation rate t_r

With t_r = 0 the whorls stay in a single plane (a planispiral); increasing it stacks the whorls into a high-spired form.

compare_raup([0.0, 1.0, 4.0], vary="t_r", title="Translation rate",
             w_r=2.0, d_r=0.0).show()

Relative distance from the coiling axis d_r

d_r (in (-1, 1)) sets how far the generating curve sits from the coiling axis.

compare_raup([-0.4, 0.0, 0.4], vary="d_r", title="Distance from the axis",
             w_r=2.0, t_r=1.0).show()

Generate shells in batch

RaupModel provides the same model in a scikit-learn-style estimator. Its inverse_transform is the generative map from parameters to form.

model = RaupModel()
params = np.array(
    [
        [1.2, 0.5, 0.1],
        [1.5, 1.5, 0.2],
        [2.0, 0.0, 0.05],
    ]
)
surfaces = model.inverse_transform(params, theta_range=theta, phi_range=phi)
surfaces.shape

(3, 240, 60, 3)

The leading axis indexes the three parameter rows; the rest is the familiar (n_theta, n_phi, 3) surface.

plot_surface3d(surfaces[1], title="Second parameter set").show()

For analysis pipelines, pass as_frame=True to get a tidy long-format DataFrame indexed by (specimen_id, trajectory_id, phi_id).

df = model.inverse_transform(
    params, theta_range=theta, phi_range=phi, as_frame=True
)
df.head()

			x	y	z
specimen_id	trajectory_id	phi_id
0	0	0	2.222222	0.0	1.111111
		1	2.216557	0.0	1.217405
		2	2.199626	0.0	1.322494
		3	2.171620	0.0	1.425188
		4	2.132857	0.0	1.524323

Sample evenly along the shell

By default the surface is sampled at equal steps in the coiling angle theta. Because the shell is a logarithmic spiral, equal theta steps pack many points near the tight apex and stretch them apart at the wide aperture. For even meshes, sample at equal steps of arc length along the shell instead.

l_r and theta_r convert between the coiling angle and the trajectory arc length l. To space n_theta samples evenly in arc length, take a uniform grid in l and map it back to theta with theta_r:

w_r, t_r, d_r = 2.0, 0.7, 0.05
n_whorls = 4.0
n_theta = 48
phi = np.linspace(0.0, 2.0 * np.pi, 40)

theta_uniform = np.linspace(0.0, 2.0 * np.pi * n_whorls, n_theta)

l_max = l_r(2.0 * np.pi * n_whorls, w_r, t_r, d_r)
theta_arclength = theta_r(np.linspace(0.0, l_max, n_theta), w_r, t_r, d_r)

fig = make_subplots(
    rows=1,
    cols=2,
    specs=[[{"type": "surface"}, {"type": "surface"}]],
    subplot_titles=["equal steps in theta", "equal steps in arc length"],
    horizontal_spacing=0.01,
)
for col, theta in enumerate([theta_uniform, theta_arclength], start=1):
    Xi = raup(w_r, t_r, d_r, theta_range=theta, phi_range=phi)
    Xi = Xi * np.array([-1.0, 1.0, -1.0])  # apex-up (see plot_surface3d note)
    fig.add_trace(
        go.Surface(
            x=Xi[..., 0],
            y=Xi[..., 1],
            z=Xi[..., 2],
            showscale=False,
            colorscale="Viridis",
        ),
        row=1,
        col=col,
    )
fig.update_layout(
    width=900,
    height=450,
    margin=dict(l=0, r=0, t=40, b=0),
    scene=dict(aspectmode="data"),
    scene2=dict(aspectmode="data"),
)
fig.show()

Notice how the mesh on the right is spaced much more evenly along the coil.

See also

• Growing tube model to generate shells with the growing tube model, including heteromorph (non-constant) growth.

• Theoretical morphological models of coiling for the theory behind Raup’s model.

References

• Raup, D.M., Michelson, A., 1965. Theoretical Morphology of the Coiled Shell. Science 147, 1294–1295.

• Raup, D.M., 1966. Geometric analysis of shell coiling: general problems. Journal of Paleontology 40, 1178–1190.