bruno/catastrophe
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

use xz plane instead of xy

Browse source
This commit is contained in:
Bruno Postle 2026-03-25 21:48:29 +00:00
parent af1c43d41c
commit 9b33b566dd
3 changed files with 316920 additions and 892709 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

47
CLAUDE.md
View file

@ -19,21 +19,22 @@ Each surface is the **equilibrium manifold** — the set of all points where the
> dV/dx = 0 > dV/dx = 0
The 3D surface is swept over the control parameter space (a, b), with x (the state variable) as the third axis. Where the surface folds back on itself is the **bifurcation set** — the region where the system can catastrophically jump between states.
### Butterfly Catastrophe ### Butterfly Catastrophe
- **Potential:** V(x) = x⁶ + ax⁴ + cx² + dx (a = −3 fixed, b = 0) - **Potential:** V(x) = x⁶ + ax⁴ + cx² + dx (a = −3 fixed, b = 0)
- **Equilibrium condition:** dV/dx = 6x⁵ − 12x³ + 2cx + d = 0 - **Equilibrium condition:** dV/dx = 6x⁵ − 12x³ + 2cx + d = 0
- **Control space:** (c, d) swept over a 2D grid - **Rearranged as a height field:** d = −6x⁵ + 12x³ − 2cx
- **State space:** up to 5 real roots x at any given (c, d) - **Print base:** (x, c) plane — state variable × control parameter
- **Characteristic feature:** nested "butterfly wing" fold structure — a self-intersecting bifurcation curve in the (c, d) plane enclosing a 5-root "pocket" (c∈[0,3], d≈0), surrounded by a 3-root wing region, with a single-root region outside - **Print height:** d (the other control parameter, computed directly)
> **Why not vary (a, b) with c=d=0?** With c=d=0, the equilibrium equation factors as This parameterisation is **single-valued**: every (x, c) point maps to exactly one d, so the mesh is a simple height field with no multi-valued branches, no root finding, and no fold-edge gaps.
> x²(6x³ + 4ax + 3b) = 0 — x=0 is always a double root and the remaining roots come
> from a cubic, which is structurally identical to the **cusp** catastrophe. The butterfly The **fold ridges** — where the surface has zero gradient in x — satisfy ∂d/∂x = 0, giving the bifurcation curve c = 18x² − 15x⁴. This self-intersecting curve is visible as a characteristic ridge on the surface.
> structure only appears when d ≠ 0 generically, which requires d (or an equivalent odd
> perturbation) to be varied as a control parameter. > **Why not sweep (c, d) and solve for x?** That approach requires finding multiple roots
> of a degree-5 polynomial at each grid point, tracking which root belongs to which branch
> across fold lines, and capping fold edges — all of which introduce artefacts and holes.
> Rearranging to d(x, c) avoids all of this entirely.
--- ---
@ -49,11 +50,10 @@ The 3D surface is swept over the control parameter space (a, b), with x (the sta
## How the Generator Works ## How the Generator Works
1. **Root finding** — at each (a, b) grid point, all real roots of dV/dx = 0 are found using Newton-Raphson with dense initial seeding across the x range 1. **Height field** — for each (x, c) grid point, compute d = (−6x⁵ + 12x³ − 2cx) × D_SCALE
2. **Branch tracking** — roots are sorted and matched by branch index across adjacent grid cells 2. **Mesh construction** — adjacent grid quads are triangulated into a regular height-field mesh
3. **Mesh construction** — adjacent grid quads on the same branch are triangulated into a surface mesh 3. **Base slab** — a flat rectangular base is added so the model is self-supporting on a print bed
4. **Base slab** — a flat rectangular base is added so the model is self-supporting on a print bed 4. **ASCII STL output** — written as ASCII (not binary) for maximum compatibility with slicers and viewers
5. **ASCII STL output** — written as ASCII (not binary) for maximum compatibility with slicers and viewers
--- ---
@ -69,13 +69,12 @@ Output: `butterfly_catastrophe.stl`
| Parameter | Default | Effect | | Parameter | Default | Effect |
|-----------|---------|--------| |-----------|---------|--------|
| `GRID` | `50` | Resolution of the (c,d) control grid — increase to 80–100 for final print | | `GRID` | `200` | Grid resolution — higher = smoother fold ridges |
| `C_RANGE` | `(-2.0, 7.0)` | Range of control parameter c | | `X_RANGE` | `(-3.0, 3.0)` | Range of state variable x (print width) |
| `D_RANGE` | `(-6.0, 6.0)` | Range of control parameter d | | `C_RANGE` | `(-1.0, 5.0)` | Range of control parameter c (print depth) |
| `X_RANGE` | `2.5` | Search window for equilibrium roots | | `D_SCALE` | `0.5` | Vertical scale factor — reduce if the model is too tall |
| `MAX_MATCH_DZ` | `0.8` | Max z-gap for inter-row branch matching in main surface |
For a final high-quality print, increase `grid` to `80`–`100`. The default of `40` is optimised for STL viewer compatibility. The butterfly fold structure is concentrated around x ∈ [−1.1, 1.1] and c ∈ [0, 5.4]; extending X_RANGE beyond ±2 adds flat outer wings with no additional features.
--- ---
@ -83,8 +82,8 @@ For a final high-quality print, increase `grid` to `80`–`100`. The default of
- **Orientation:** flat base down — no supports needed - **Orientation:** flat base down — no supports needed
- **Layer height:** 0.15–0.20 mm for good surface detail - **Layer height:** 0.15–0.20 mm for good surface detail
- **Perimeters:** ≥ 2, as the fold regions are thin - **Perimeters:** ≥ 2 for the thin ridge regions
- **Scale:** ~120 mm along the a-axis makes a good desk model - **Scale:** the fold ridges are most visible at ~100–150 mm along the c-axis
- **Material:** PLA or PETG both work well; the overhangs are gentle - **Material:** PLA or PETG both work well; the overhangs are gentle
--- ---
@ -95,4 +94,4 @@ For a final high-quality print, increase `grid` to `80`–`100`. The default of
numpy numpy
``` ```
No other dependencies — STL writing uses Python's built-in `struct` module (binary) or plain file I/O (ASCII). No other dependencies — STL writing uses plain file I/O (ASCII STL).

298
butterfly_catastrophe.py
View file

@ -4,285 +4,93 @@ Butterfly Catastrophe Surface — ASCII STL Generator
Potential: V(x) = x^6 + a*x^4 + c*x^2 + d*x (a = -3 fixed, b = 0) Potential: V(x) = x^6 + a*x^4 + c*x^2 + d*x (a = -3 fixed, b = 0)
Equilibrium: dV/dx = 6x^5 - 12x^3 + 2c*x + d = 0 Equilibrium: dV/dx = 6x^5 - 12x^3 + 2c*x + d = 0
Control parameters (print base): c (horizontal), d (depth) Rearranged as a single-valued function:
State variable (print height): x d(x, c) = -(6x^5 - 12x^3 + 2c*x)
= -6x^5 + 12x^3 - 2c*x
With a = -3, the bifurcation set in the (c, d) plane forms the characteristic The print base is the (x, c) plane; height is d. Every grid point maps to
butterfly shape: a self-intersecting loop passing through (c=3, d=0), enclosing exactly one vertex — no root finding, no branch tracking, no holes.
a region with 5 equilibria ("butterfly pocket"), surrounded by a 3-root region
with outer fold wings, and a single-root region outside. This is structurally The fold lines of the bifurcation set appear as ridges where the surface has
different from the cusp catastrophe and cannot be seen with c = d = 0 fixed. zero gradient in x: ∂d/∂x = -30x^4 + 36x^2 - 2c = 0, i.e. c = 18x^2 - 15x^4.
""" """
import numpy as np import numpy as np
import os import os
A_FIXED = -3.0 # butterfly unfolding parameter (must be negative) A_FIXED = -3.0
# ── Tuning ────────────────────────────────────────────────────────────────── # ── Tuning ──────────────────────────────────────────────────────────────────
GRID = 200 # control-space resolution — increase to 100–120 for print GRID = 200 # grid resolution — higher = smoother ridges
C_RANGE = (-2.0, 7.0) # range of control parameter c X_RANGE = (-3.0, 3.0) # state variable x
D_RANGE = (-6.0, 6.0) # range of control parameter d C_RANGE = (-1.0, 5.0) # control parameter c
X_RANGE = 2.5 # half-width of root search window D_SCALE = 0.5 # scale factor applied to the computed d height;
MAX_MATCH_DZ = 0.8 # max z-gap for inter-row branch matching # reduce if the model is too tall for your printer
# ── 1. Root finding ────────────────────────────────────────────────────────── # ── 1. Height function ────────────────────────────────────────────────────────
def dV(x, c, d): def d_surface(x, c):
return 6*x**5 + 4*A_FIXED*x**3 + 2*c*x + d """d value on the equilibrium manifold: dV/dx = 0 solved for d."""
return (-6*x**5 + 12*x**3 - 2*c*x) * D_SCALE
def d2V(x, c, d):
return 30*x**4 + 12*A_FIXED*x**2 + 2*c
def find_roots(c, d, n_starts=80):
"""Return sorted real roots of dV/dx = 0 for the given (c, d)."""
xs = np.linspace(-X_RANGE, X_RANGE, n_starts)
roots = []
for x0 in xs:
x = float(x0)
for _ in range(200):
fx = dV(x, c, d)
if abs(fx) < 1e-12:
break
dfx = d2V(x, c, d)
if abs(dfx) < 1e-14:
break
step = fx / dfx
x -= step
if abs(x) > 2.0 * X_RANGE: # diverged — abandon
break
if abs(step) < 1e-10:
break
if abs(dV(x, c, d)) < 1e-7 and abs(x) <= X_RANGE + 0.15:
if not any(abs(x - r) < 1e-4 for r in roots):
roots.append(x)
return sorted(roots)
# ── 2. Branch tracking ───────────────────────────────────────────────────────
def track_branches(roots_along_axis):
"""
Track branches along one axis (fixed d, varying c) using greedy
nearest-neighbour matching. Returns a list of tracks; each track is a
list of length GRID where entry i is the root value at column i, or None
when the branch does not exist there.
Sorting-index matching (branch 0 always = branch 0) breaks at folds
because two adjacent branches coalesce, shifting every higher index by
one. Nearest-neighbour tracking follows each physical sheet through the
fold correctly: the two merging branches each get None past the fold, and
the surviving sheet keeps its track unbroken.
"""
n = len(roots_along_axis)
if n == 0:
return []
tracks = [[r] for r in roots_along_axis[0]]
for i in range(1, n):
curr = roots_along_axis[i]
prev_live = [(ti, t[-1]) for ti, t in enumerate(tracks)
if t[-1] is not None]
prev_matched, curr_matched = set(), set()
matches = {} # track_idx → curr_root_idx
cands = sorted(
[(abs(pv - curr[ci]), ti, ci)
for ti, pv in prev_live
for ci in range(len(curr))],
key=lambda x: x[0]
)
for _, ti, ci in cands:
if ti not in prev_matched and ci not in curr_matched:
matches[ti] = ci
prev_matched.add(ti)
curr_matched.add(ci)
for ti, t in enumerate(tracks):
t.append(curr[matches[ti]] if ti in matches else None)
# Branches that appear for the first time at this column
for ci in range(len(curr)):
if ci not in curr_matched:
tracks.append([None] * i + [curr[ci]])
return tracks
# ── 3. Build mesh ────────────────────────────────────────────────────────────
def _emit_quad(triangles, p00, p10, p11, p01):
triangles.append((p00, p10, p11))
triangles.append((p00, p11, p01))
def _fold_terminations(tracks, axis_vals):
"""
Scan tracks along one axis and return a list of fold-termination events.
Each event is (axis_val, xa, xb) where axis_val is the last valid position,
and xa < xb are the two branch values that die together at a fold.
Branches come in pairs at fold lines (two coalesce), so we pair adjacent
sorted dying values. Events are indexed by the axis position of the last
valid step.
"""
events = []
n = len(axis_vals) - 1 # number of steps
for step in range(n):
dying = sorted(t[step] for t in tracks
if t[step] is not None and t[step + 1] is None)
for k in range(0, len(dying) - 1, 2):
events.append((axis_vals[step], dying[k], dying[k + 1]))
return events
# ── 2. Build mesh ─────────────────────────────────────────────────────────────
def build_mesh(): def build_mesh():
x_vals = np.linspace(*X_RANGE, GRID)
c_vals = np.linspace(*C_RANGE, GRID) c_vals = np.linspace(*C_RANGE, GRID)
d_vals = np.linspace(*D_RANGE, GRID)
print(f' Computing roots on {GRID}×{GRID} grid…') # Pre-compute the full height field in one vectorised call
# roots_grid[i][j] = sorted roots at (c_vals[i], d_vals[j]) X, C = np.meshgrid(x_vals, c_vals, indexing='ij') # (GRID, GRID)
roots_grid = [[find_roots(c, d) for d in d_vals] for c in c_vals] D = (-6*X**5 + 12*X**3 - 2*C*X) * D_SCALE
# Track branches along rows (fixed d, varying c) and columns (fixed c, varying d).
print(' Tracking branches…')
row_tracks = [track_branches([roots_grid[i][j] for i in range(GRID)])
for j in range(GRID)]
col_tracks = [track_branches([roots_grid[i][j] for j in range(GRID)])
for i in range(GRID)]
triangles = [] triangles = []
# ── Main surface quads ────────────────────────────────────────────────────
for j in range(GRID - 1):
d0, d1 = d_vals[j], d_vals[j + 1]
tracks_j = row_tracks[j]
tracks_j1 = row_tracks[j + 1]
for i in range(GRID - 1):
c0, c1 = c_vals[i], c_vals[i + 1]
segs_j = [(t[i], t[i + 1], k) for k, t in enumerate(tracks_j)
if t[i] is not None and t[i + 1] is not None]
segs_j1 = [(t[i], t[i + 1], k) for k, t in enumerate(tracks_j1)
if t[i] is not None and t[i + 1] is not None]
if not segs_j or not segs_j1:
continue
j1_used = set()
for x00, x10, _ in sorted(segs_j, key=lambda s: (s[0] + s[1]) / 2):
best_k1 = best_x01 = best_x11 = None
best_dist = MAX_MATCH_DZ
for x01, x11, k1 in segs_j1:
if k1 in j1_used:
continue
dist = max(abs(x00 - x01), abs(x10 - x11))
if dist < best_dist:
best_dist, best_k1 = dist, k1
best_x01, best_x11 = x01, x11
if best_k1 is None:
continue
j1_used.add(best_k1)
_emit_quad(triangles,
(c0, d0, x00), (c1, d0, x10),
(c1, d1, best_x11), (c0, d1, best_x01))
# ── C-direction fold caps ─────────────────────────────────────────────────
# The fold line runs at an angle through the (c, d) grid, so the column
# where two branches die can differ by ±1 between adjacent rows.
# We collect all fold terminations per row, then match them across the
# row pair regardless of exact column, connecting the dying edges with
# (possibly trapezoidal) quads.
c_terms = [_fold_terminations(row_tracks[j], c_vals) for j in range(GRID)]
for j in range(GRID - 1):
d0, d1 = d_vals[j], d_vals[j + 1]
terms_j = c_terms[j]
terms_j1 = c_terms[j + 1]
j1_used = set()
for c_j, xa, xb in terms_j:
mid = (xa + xb) / 2
best_k = None
best_dist = 1.0 # max allowed x-midpoint distance between matched pairs
for k1, (c_j1, xa1, xb1) in enumerate(terms_j1):
if k1 in j1_used:
continue
dist = abs(mid - (xa1 + xb1) / 2)
if dist < best_dist:
best_dist, best_k = dist, k1
if best_k is None:
continue
j1_used.add(best_k)
c_j1, xa1, xb1 = terms_j1[best_k]
# Cap quad: lies at the fold edge, spanning d0→d1 between the two
# dying branches. c may differ slightly between the two rows if
# the fold line is diagonal.
_emit_quad(triangles,
(c_j, d0, xa), (c_j, d0, xb),
(c_j1, d1, xb1), (c_j1, d1, xa1))
# ── D-direction fold caps ─────────────────────────────────────────────────
d_terms = [_fold_terminations(col_tracks[i], d_vals) for i in range(GRID)]
for i in range(GRID - 1): for i in range(GRID - 1):
c0, c1 = c_vals[i], c_vals[i + 1] for j in range(GRID - 1):
terms_i = d_terms[i] p00 = (x_vals[i], c_vals[j], D[i, j ])
terms_i1 = d_terms[i + 1] p10 = (x_vals[i+1], c_vals[j], D[i+1, j ])
p11 = (x_vals[i+1], c_vals[j+1], D[i+1, j+1])
i1_used = set() p01 = (x_vals[i], c_vals[j+1], D[i, j+1])
for d_i, xa, xb in terms_i: triangles.append((p00, p10, p11))
mid = (xa + xb) / 2 triangles.append((p00, p11, p01))
best_k = None
best_dist = 1.0
for k1, (d_i1, xa1, xb1) in enumerate(terms_i1):
if k1 in i1_used:
continue
dist = abs(mid - (xa1 + xb1) / 2)
if dist < best_dist:
best_dist, best_k = dist, k1
if best_k is None:
continue
i1_used.add(best_k)
d_i1, xa1, xb1 = terms_i1[best_k]
_emit_quad(triangles,
(c0, d_i, xa), (c1, d_i1, xa1),
(c1, d_i1, xb1), (c0, d_i, xb))
return triangles return triangles
# ── 4. Flat base ───────────────────────────────────────────────────────────── # ── 3. Flat base ──────────────────────────────────────────────────────────────
def add_base(triangles, z_base=-2.8): def add_base(triangles):
x0, x1 = X_RANGE
c0, c1 = C_RANGE c0, c1 = C_RANGE
d0, d1 = D_RANGE X, C = np.meshgrid(np.linspace(x0, x1, GRID),
zb, zt = z_base, z_base + 0.15 np.linspace(c0, c1, GRID), indexing='ij')
D = (-6*X**5 + 12*X**3 - 2*C*X) * D_SCALE
zb = D.min() - 0.15
zt = zb + 0.15
# Top and bottom faces of base slab
triangles += [ triangles += [
((c0,d0,zt),(c1,d0,zt),(c1,d1,zt)), ((x0,c0,zt),(x1,c0,zt),(x1,c1,zt)),
((c0,d0,zt),(c1,d1,zt),(c0,d1,zt)), ((x0,c0,zt),(x1,c1,zt),(x0,c1,zt)),
((c0,d0,zb),(c1,d1,zb),(c1,d0,zb)), ((x0,c0,zb),(x1,c1,zb),(x1,c0,zb)),
((c0,d0,zb),(c0,d1,zb),(c1,d1,zb)), ((x0,c0,zb),(x0,c1,zb),(x1,c1,zb)),
] ]
for (x0,y0),(x1,y1) in [ for (ax,ay),(bx,by) in [
((c0,d0),(c1,d0)), ((c1,d0),(c1,d1)), ((x0,c0),(x1,c0)), ((x1,c0),(x1,c1)),
((c1,d1),(c0,d1)), ((c0,d1),(c0,d0)), ((x1,c1),(x0,c1)), ((x0,c1),(x0,c0)),
]: ]:
triangles += [ triangles += [
((x0,y0,zb),(x1,y1,zb),(x1,y1,zt)), ((ax,ay,zb),(bx,by,zb),(bx,by,zt)),
((x0,y0,zb),(x1,y1,zt),(x0,y0,zt)), ((ax,ay,zb),(bx,by,zt),(ax,ay,zt)),
] ]
return triangles return triangles
# ── 5. ASCII STL output ─────────────────────────────────────────────────────── # ── 4. ASCII STL output ───────────────────────────────────────────────────────
def normal(v0, v1, v2): def normal(v0, v1, v2):
a = np.subtract(v1, v0) a = np.subtract(v1, v0)
b = np.subtract(v2, v0) b = np.subtract(v2, v0)
n = np.cross(a, b) n = np.cross(a, b)
length = np.linalg.norm(n) L = np.linalg.norm(n)
return n / length if length > 1e-14 else np.array([0.0, 0.0, 1.0]) return n / L if L > 1e-14 else np.array([0.0, 0.0, 1.0])
def write_ascii_stl(triangles, filename): def write_ascii_stl(triangles, filename):
with open(filename, 'w') as f: with open(filename, 'w') as f:
@ -301,10 +109,10 @@ def write_ascii_stl(triangles, filename):
print(f' Triangles : {len(triangles):,}') print(f' Triangles : {len(triangles):,}')
print(f' File size : {size_kb:.0f} KB → {filename}') print(f' File size : {size_kb:.0f} KB → {filename}')
# ── 6. Main ─────────────────────────────────────────────────────────────────── # ── 5. Main ───────────────────────────────────────────────────────────────────
if __name__ == '__main__': if __name__ == '__main__':
print(f'Building butterfly catastrophe mesh (grid={GRID})…') print(f'Building butterfly catastrophe surface (grid={GRID})…')
tris = build_mesh() tris = build_mesh()
tris = add_base(tris) tris = add_base(tris)
print('Writing ASCII STL…') print('Writing ASCII STL…')

1209284
butterfly_catastrophe.stl
View file

File diff suppressed because it is too large Load diff