Mandelbulb

DEFINITION

vₖ₊₁ = vₖⁿ + c,   v₀ = 0,   c ∈ ℝ³
vⁿ = rⁿ(cos nθ·cos nφ,  sin nθ·cos nφ,  −sin nφ),   r = |v|
θ = atan2(y, x),   φ = atan(z / √(x²+y²)),   n = 8
dr → n·rⁿ⁻¹·dr + 1,   d = r·ln r / (2·dr)   — heuristic, not a bound
sphere trace: t → t + 0.75·d(o + t·û),   n̂ = ∇d / |∇d|

NOTES

There is no three-dimensional number system to run the Mandelbrot recursion in. Hamilton spent years failing to build one and reached the quaternions in four dimensions instead, and that no such three-dimensional algebra exists is a theorem, so every honest attempt at a solid Mandelbrot has had to invent something that is not one. What Daniel White proposed on a fractal forum in November 2007, and Paul Nylander generalised to arbitrary powers in August 2009, is the cheerful workaround: forget multiplication and copy what the complex square DOES to a point in polar form, which is to square the length and double the angle. Read a point in spherical coordinates, raise the radius to the nth power, multiply BOTH angles by n, and add the constant. It is not associative, it is not distributive, it is barely a product at all, and iterating it from zero gives a solid with structure in all three directions, which is what the thread it was posted in had spent two years failing to get out of quaternions: the complaint there, in its own words, is that a quaternion fractal has Mandelbrot-depth complexity in only one dimension and looks smooth on the other two. The power that made it famous is eight. What this plate draws is that solid, lit, found by marching rays at it rather than by cutting a plane through it. The march needs a distance, and where the distance comes from is where the honesty of this plate lives. Carry one extra number beside the orbit, sent at every step to n times the radius to the n minus one times itself plus one, and the radius times its own logarithm over twice that number is the estimate of the distance to the boundary that descends from the Douady-Hubbard potential of a filled Julia set. That is a theorem for a map that is holomorphic, and this one is not. The triplex power multiplies lengths exactly — measured here to fifteen figures, and it is exactly what carries the escape argument across from the plane — but a length is not a Jacobian. The map multiplies both spherical angles by n, so near the poles the azimuthal circle has a radius falling to zero while its image is traversed n times, and the tangential stretch there grows without bound where the radial factor does not. The single number therefore understates how fast the map can stretch, so the estimate overstates how far a ray may safely travel, so a march can step through the surface. It is a heuristic and not a bound, and the guarantee that makes sphere tracing exact does not apply to it. The mitigation is the one the published lineage gives and it is a trade rather than a cure: take three quarters of the offered step. What that buys was measured rather than assumed, and the answer is not the obvious one. The silhouette is stable at every relaxation factor tested, down to a handful of pixels in a frame; what improves is where along the ray the march stops, and that is what the reader is actually looking at when the valley floors between the equatorial lobes come out as soft washes where the structure should keep subdividing, and when fine concentric rings appear in the ribbed whorls. Both are the estimate stopping short of the surface. Both are still there at the factor that ships, and the numbers are under Provenance rather than smoothed over. Where the mathematics stops and the picture begins has to be said plainly, because a lit render invites the confusion. The estimate, the march, where it stops and the normal taken from the gradient of the same estimate are the plate. The light is not: two directions fixed relative to the eye with one slider turning the key around, and where a lamp stands is staging rather than a fact about the set. Lambert on that normal, a halfway-vector specular, a cone traced back toward the light for the soft shadow and a hemisphere of short rays for the ambient obscurance are published models rather than invented ones, which is why they are named under Provenance and not merely used. And the optics of the page itself run after all of it and outside this function entirely, exactly as they do over every other plate here. The camera does not move, and it cannot. Rendering this way costs some tens of times a plane cut, far more than one frame holds, so the picture develops instead: every paint traces a batch of rays into a buffer, coarsest first, the first paint laying down a whole eighth-resolution picture and later ones refining it until every cell has been traced once. A camera move would throw all of that away, so the view is two constants the reader turns rather than an animation, and turning one restarts the development, which takes about a second. The shader does the same march in one frame at full resolution and arrives at the same picture, and that it is the same picture is measured under Provenance rather than assumed. The frame is set by the object rather than chosen: the containment radius two to the power one over n minus one bounds this solid for the same reason it bounds the complex multibrot, since the triplex power multiplies lengths exactly, and the camera sits at 2.9 of it, so every power frames itself. That radius is why the power slider carries the history of the plate without any other dial moving. Set it to eight for the object White and Nylander rendered and the name attached to; set it to two and it is the triplex squaring White posted in 2007, whose mass is concentrated well inside its own bound while thin filaments reach almost to it, so that object sits smaller in the frame than the others — a fact about it, measured, not a framing error. And the identity the previous edition of this plate was built on is still true and is no longer what the picture shows. Hold the third coordinate at zero and it stays zero, exactly: the elevation angle is zero, its cosine is one and its sine is zero, so the equatorial plane is invariant and the map restricted to it is the complex nth power written out. The equatorial section of this solid is therefore the degree-n multibrot, computed by an identical recursion, and at n = 2 it is the Mandelbrot set itself. That was checked at the level of the orbit and it stands. What cannot be claimed is that the picture on this plate shows it: a lit solid seen from outside shows a surface, and the section that carries the identity is a plane through the middle of it that no ray reaches. The claim has moved from the picture to the mathematics, which is where it was always true. Rudy Rucker had tried nearly this squaring in spherical coordinates in 1988 and set it aside, and came back to write about it in 2009 while the forum thread was running. None of the three ever published this construction in a journal.

PROVENANCE

Origin

Daniel White, posting as "twinbee", fractalforums.com thread "True 3D mandlebrot type fractal", opened 22 September 2007. The triplex squaring is his post #10, 17 November 2007: r = sqrt(x*x+y*y+z*z), yang = atan2(sqrt(x*x+y*y), z), zang = atan2(y, x), then newx = r*r*sin(yang*2 + 0.5*pi)*cos(zang*2 + pi), newy the same with sin(zang*2 + pi), newz = r*r*cos(yang*2 + 0.5*pi). Read off the archived post itself, not from a secondary account. The misspelling in the thread title is the original: the topic line was edited on 19 October 2009 and every reply subject still carries "mandlebrot".

The power-n form this plate implements

Paul Nylander, posting as "bugman", same thread, post #153, 30 August 2009, quoted verbatim from the archive: "{x,y,z}^n = r^n{cos(n*theta)cos(n*phi),sin(n*theta)cos(n*phi),-sin(n*phi)}, r=sqrt(x^2+y^2+z^2), theta=atan(y/x), phi=atan(z/sqrt(x^2+y^2))". Two conventions differ here and the difference is not cosmetic: White measures his angle from the vertical axis and carries two phase offsets, Nylander measures it from the horizontal plane and carries none. Worked through, the two offsets White carries are exactly what converts one convention into the other at n = 2, so the two agree there and only there; carried naively to n = 8 they give the map of Nylander composed with a half turn. This plate implements the form Nylander posted, at every power, so power 2 reproduces the object White posted exactly rather than up to a symmetry. The angle Nylander writes as atan(z/sqrt(x^2+y^2)) is the same angle as asin(z/r), since r squared is rho squared plus z squared with rho never negative, and the implementation calls neither function: it uses the cosine and sine of that angle directly, as the pair (sqrt(x^2+y^2)/r, z/r). The vertical axis of the map is the vertical axis of the render, so the poles the theory warns about are at the top and bottom of the frame and the reader can go and look at them.

Never published by its discoverers

No journal paper by White or Nylander exists for this object. It was published in the plain sense of being posted in the open, in a public thread anyone could read and argue with, and that thread is what src points at. The archive at ff1.fractalforums.org preserves the original topic numbering; the old host www.fractalforums.com answers with HTTP 200 and a database error in the body, and is not linked from here for that reason.

The one it is not

A claim that Jules Ruis constructed this object in 1997 circulates online and is repeated in a number of secondary write-ups. It does not stand up to the checking this atlas does before it prints a name. No 1997 document showing the formula resolves; the earliest reachable material making the claim is a 2018 piece by Ruis himself, which describes a 1997 BASIC program called BBM-15.exe rather than reproducing any construction. None of the people involved credits him: not the article by White, not the post by Nylander, not the blog post by Rucker, and not the one academic treatment of the object, which mentions him nowhere. English Wikipedia does not carry the claim either, crediting White and Nylander in 2009. This paragraph exists because a wrong attribution here would be worse than none, and because a reader who has met the claim elsewhere deserves to know it was looked at.

The precursor, and the name

Rudy Rucker had tried a spherical-coordinate squaring around 1988 and set it aside, and returned to it in "In Search of the Mandelbulb" on his own blog, 15 September 2009, where he writes plainly: "This is an idea I spent a little time on in 1988, and recently I came back to it again." Who coined the name is a small muddle nobody needs resolved. Rucker credits White in that same post — "And White coined the world [sic] Mandelbulb for the image we seek" — while Nylander credits Rucker. The honest summary is that the name emerged between the three of them in the autumn of 2009; it is in casual use by strangers in the forum thread itself by 11 November 2009, alongside the article by White, "The Unravelling of the Real 3D Mandelbrot Fractal", 8 November 2009.

Post-hoc academic treatment

V. Boily and D. Rochon, "On the algebraic foundation of the Mandelbulb", Fractals 31(05), 2023, doi:10.1142/S0218348X23500627 — verified against the live Crossref record on title, both authors, journal, volume, issue and year, and read in full. They build a spherical product on the pure quaternions, show the structure is a commutative unital magma, and present a set they describe as visually identical to the Mandelbulb, proving it bounded and generable by an escape-time algorithm with the same boundary radius the multibrots have. One thing this plate deliberately does NOT take from them: their abstract says one of its two-dimensional cuts has the same dynamics as the Mandelbrot set, and that sentence is about a different construction in the same paper, their bulbic Mandelbrot set, not about the Mandelbulb. The equatorial-section claim in the note above is this plate own, and was measured here.

The rendering, and every part of it named

The march is sphere tracing: J. C. Hart, "Sphere tracing: a geometric method for the antialiased ray tracing of implicit surfaces", The Visual Computer 12(10), 1996, 527-545, doi:10.1007/s003710050084. Given a function returning a lower bound on the distance to the surface, step by exactly that bound; each evaluation defines an unbounding sphere guaranteed to contain no intersection, so no step can penetrate. What this plate has is not that bound, which is the subject of the next block. The soft shadow is section 5 of the same paper, not a folk trick: marching toward the light, the unbounding sphere of radius d at distance t subtends a half-angle whose sine is d/t, so the running minimum of that ratio bounds how much of the light cone is still visible. The estimate itself is the Douady-Hubbard form, stated in Peitgen and Saupe (eds), The Science of Fractal Images, Springer 1988, doi:10.1007/978-1-4612-3784-6, and descending from Douady and Hubbard, Etude dynamique des polynomes complexes (Orsay notes 1984-85, genuinely pre-DOI and rendering no link). The lineage of distance-bounded marching on a deterministic fractal is J. C. Hart, D. J. Sandin and L. H. Kauffman, "Ray tracing deterministic 3-D fractals", ACM SIGGRAPH Computer Graphics 23(3), 1989, 289-296, doi:10.1145/74334.74363. Shading: Lambert, Photometria, 1760, pre-DOI and public domain, of which the whole content here is max(0, n dot l); the specular is the halfway vector of J. F. Blinn, "Models of light reflection for computer synthesized pictures", Proceedings of SIGGRAPH 77, 192-198, doi:10.1145/563858.563893 (the exponent-on-the-reflection alternative is Bui Tuong Phong, "Illumination for computer generated pictures", Communications of the ACM 18(6), 1975, 311-317, doi:10.1145/360825.360839, and is one line either way). The ambient term is S. Zhukov, A. Iones and G. Kronin, "An ambient light illumination model", Eurographics Rendering Techniques 98, Springer 1998, 45-55, doi:10.1007/978-3-7091-6453-2_5: W(P) is the mean over the hemisphere of a monotone function of the distance to the nearest occluder, zero at zero and one at infinity, and what runs here is four cosine-weighted directions each sphere traced to a radius R with the simplest member of the family they give. It is deliberately NOT the five-tap distance-field trick that circulates as ambient occlusion in demo code: that has no published derivation, and reproducing it would be reproducing a sketch. These are the same six identifiers mandelbox (PL. 149) carries, all resolved against the live Crossref API on 26 August 2026 and matched on title, authors, container, year and pages; Blinn 1977 carries empty volume and issue fields there, which is why it is cited as the proceedings rather than in the Computer Graphics 11(2) form a secondary source would supply.

The estimate is a heuristic, and here is what that costs

mandelbox can prove its estimate is a Lipschitz bound, because every factor of that map is a similarity. This one cannot: the triplex power multiplies lengths exactly but multiplies both spherical angles by n, so the tangential stretch is unbounded near the poles while the scalar tracks only the radial factor. What was measured, rather than argued. Penetration, the same test mandelbox ran on itself: every full march step subdivided 24 ways and every sub-point tested for membership of the same 12-iteration bounded set the estimate is estimating, at four views spanning both the equator and the pole and at powers 2, 8 and 12. At the relaxation factor that ships, 0 of 69,031 steps penetrate at n = 8 near the equator, 1 of 60,860 at n = 8 looking down the pole (and that one at 96 per cent of the way through its own step), 0 of 57,324 at n = 12, 0 of 49,875 at n = 12 down the pole, and 5 of 64,584 at n = 2, which is the worst case at 0.008 per cent. Just under 7.2 million membership tests in total. So the theoretical failure is real and is rare here, and the polar view — where the mathematics predicts the worst behaviour — is not where it shows up, which is the opposite of what the theory warns and is said here rather than quietly enjoyed. Erosion, against a reference march taking quarter steps with a budget of 2,400: the silhouette differs on 0.005 per cent of the frame at n = 8, 0.006 at n = 12 and 0.121 at n = 2, and every differing pixel is a loss rather than a gain — a bite out of the object, never an addition to it. The relaxation factor barely moves that: at full steps the same figures are 0.009, 0.022 and 0.171 per cent. What the factor does move is the hit distance, which is what the reader sees as the smoothed valley floors and the concentric rings: the mean error along the ray falls from 1.40e-3 at full steps to 9.91e-4 at three quarters, 7.15e-4 at three fifths and 5.66e-4 at a half, roughly linearly and with no knee to aim at. Three quarters ships, and the reason is the step budget rather than taste: below it the march begins to exhaust its own budget at the silhouette (2 rays of 32,200 at three fifths of a step on the card grid, more at a half), and a ray that runs out of budget renders as background, which turns a soft valley into a hard hole. The frame cost of the choice was measured too and is small: a 318 by 256 card frame on the software renderer is 23.4 ms at full steps, 24.9 at 0.85, 25.9 at three quarters, 27.8 at three fifths and 30.1 at a half.

The framing, and why it is derived rather than fixed

The camera sits at 2.9 times the containment radius two to the power one over n minus one, on a 45 degree field. That radius is a real bound and not a convenience: the triplex power multiplies lengths exactly, so |c| to the n minus one greater than 2 makes the second orbit point longer than c and the orbit runs away, exactly as in the plane. It was re-measured for the solid rather than carried over from the section, by tracing 12,000 rays inward from a far sphere in Fibonacci directions with the estimate this plate marches with: at every power from 2 to 12 the surface is found in every direction but at most two, and the furthest it reaches is 0.891 of the radius at n = 2, 0.953 at n = 3 and between 0.99 and 0.999 at every power from 4 up. So one expression frames every power, and the only visible consequence is at n = 2, where the mass of the object sits well inside a bound that only its thinnest filaments approach and it therefore reads smaller in the frame. Looked at against 2.3, 2.6, 2.9 and 3.2 of that radius at two elevations: 2.6 fills the frame to its edges with no air, 3.2 leaves the object small, and 2.9 is where the silhouette reads whole with ink around it. It is also cheaper the further out it goes, because less of the frame is a hit — a card frame on the software renderer is 33.0 ms at 2.6 with 52.2 per cent of pixels lit, 25.9 at 2.9 with 41.0 per cent, 21.7 at 3.1 with 35.5 and 17.6 at 3.4 with 29.2 — so this is a quality dial that was settled by looking and then priced, not the other way round.

What one paint can afford

One paint cannot afford a frame of this, so a paint is a batch. The grid this path shades is capped at 200 cells across, which is 32,200 rays at a 318 by 256 card and 54 batches of 600. Measured in Chrome on this machine: a batch is 2.80 ms at the median and 5.90 at its worst on a cold page, 152 ms of work in total, and 0.1 ms a frame for ever afterwards because a developed picture is a blit. The drawer at 890 by 593 is 26,600 rays in 45 batches at a median of 1.80 ms, 76 ms in total; a 240-square thumbnail is 40,000 rays in 65 batches at 4.30 ms; the 2048-wide export is 45 batches. That cap is a constraint rather than a preference: renderStill develops an exposure over 241 draws, so a resolution whose full pass needed more batches than that would ship a half-traced thumbnail and a half-traced export, which is the failure this plate would be least likely to notice. The batch is not a scanline sweep. It is hierarchical, eighth-resolution blocks first and each level halving, and the first batch is always a whole eighth-resolution pass — 525 rays at a card, inside the 600 — so a reader moving a slider sees a coarse lit solid track it within one frame rather than a black rectangle growing a corner. Every ray a coarse level traces is a ray the finished picture keeps, because a cell is sampled at its own centre and the block it fills is a nearest-neighbour spread of that one sample; sampling the centre of the block instead is the defect the mandelbox pilot found in itself, and it is invisible to everything except a comparison against the shader.

Constants, measured

Re-swept from nothing, because the constants are new and three of them mean something different. power and iter keep their keys and their meaning; the three that governed the cutting plane are gone with it, so height, tilt and drift were repurposed in place as the view elevation, the view azimuth and the key light azimuth, and the march step budget was appended. params[] is a wire format and append-only once shipped, so an old bench link still restores the power and the iteration ceiling onto the same two sliders and lands three camera angles on values they never meant, which is strictly better than the silent remap that deleting three entries would have caused. The price is a param list whose order reads oddly, and this says so. power is the one dial that changes which published object is on screen and every integer from 2 to 12 renders one: 8 is the figure White and Nylander rendered, 2 is the triplex squaring White posted in 2007, and the containment radius carries the framing across all of them without a second dial. It costs what the object costs — 11.7 ms a card frame at n = 2, 22.4 at n = 5, 25.5 at n = 8 and 30.2 at n = 12 on the software renderer — because a higher power is a fuller solid and more of the frame is a hit. iter is locked from regenerate and its old justification is void: on the plane slice it was an unusually pure cost dial, and on the solid it is not a cost dial at all. Measured across its whole range at a card frame: 24.9 ms at 4, 25.8 at 8, 25.3 at 12 and 25.6 at 16 — flat inside the noise, because the eighth power escapes so violently that a point outside the set leaves the bailout radius in two or three steps whatever the ceiling says, and the march never evaluates deep inside. What it decides instead is the definition: how many steps of the recursion get to say a point is in the set, which sets how thin a filament the estimate can resolve. It settles, and the settling is the range: against a ceiling of 24 at drawer size, a ceiling of 4 differs on 6.55 per cent of pixels, 6 on 0.290, 8 on 0.052, 12 on 0.018 and 16 on 0.014. 12 is the default and 4 to 16 is the declared range, so the bottom of it is a visibly coarser solid and the top is the reader confirming that it has converged. steps is locked for the same reason mandelbox locks it and turns out to be the same non-lever: 25.5, 25.4, 25.8 and 25.7 ms at 96, 128, 160 and 192 on the card, because the march is bound by the cost of distance estimates and not by the tail of long rays. Since it is free, it is set where the picture is settled rather than where it is cheapest: at drawer size a budget of 64 differs from 192 on 2.076 per cent of pixels, 96 on 0.441, 128 on 0.126 and 160 on 0.032, all of it at the silhouette, so 160 ships. elevation, azimuth and the light azimuth are free for regenerate and are the two least fragile constants on the plate for mandelbox own reason: a camera that always looks at the object from a distance derived from the object cannot be jittered into looking at nothing.

Checked, not just plotted

What survives from the plane slice unchanged, because it is about the map rather than about the picture, and what is new for the solid. Surviving: the norm — this power multiplies lengths exactly, since its output is the nth power of the radius times a unit vector, worst relative error 5.2e-15 over 200,000 evaluations at powers 2, 3, 5, 8 and 12, which is what carries the complex escape argument across and gives the containment radius directly. The invariance of the equatorial plane at orbit level: over 250,000 sample points per power at n = 2, 3, 8 and 12, the largest third coordinate reached by any orbit point at any step was zero — not small, zero. And the identity itself at orbit level: the same 250,000 points per power compared against an independently written complex z to the n plus c iteration gave zero escape-count mismatches at every power, with the largest difference in the smooth dwell 2.3e-10 at n = 2 and 5.3e-7 at n = 8. Retired with the section: two picture-level checks that diffed the old draw() against an independent complex multibrot renderer. They were correct and they were about a plane cut this plate no longer makes, so they are recorded here as history rather than restated as current evidence — that distinction is the whole reason the note above moves the multibrot claim from the picture to the mathematics. New for the solid: the estimate and its erosion, in their own block above; the two rendering paths, in the block below; and liveness under the real jitter. That sweep is shaped differently from every other plate in this atlas and the difference is the point: a progressive plate has no frame to sample, so each render is developed to completion rather than run for a fixed number of frames. 30 REGENERATE seeds, every constant jittered as the runner jitters it, at both a card and a drawer size, each developed until the last cell is traced: 60 renders, 0 throws, 0 below the contrast floor, worst pixelRange 131.6, median 163, best 170.1 against the gate floor of 14. The narrowest of them is n = 2, which is the dimmest object in the range for the reason the framing block gives, and it clears the floor by an order of magnitude.

The two rendering paths, checked

Both paths run the same march and they are compared rather than assumed to agree, at 20 parameter tuples: every integer power from 2 to 12 at the published constants, plus nine with both view angles, the light, the iteration ceiling and the march budget drawn at random. At a matched grid — both paths at 200 by 161, cell for cell, 1,932,000 channels — the mean channel difference is 0.151 levels in 255, 1.33 per cent of channels differ by more than one, and 0.050 per cent of pixels disagree about whether they hit the surface at all. Where that residue lives was isolated rather than asserted, by a cheaper experiment than the one mandelbox ran on itself: hold the camera at angles whose sine and cosine are exact in both float widths and it collapses. At elevation 0, azimuth 0 and the key at 0 the mean difference is 0.074 with a maximum of 42, 0.01 per cent of channels over one level and no silhouette disagreement at all; at elevation 0, azimuth 90 and the key at 180 it is 0.028 with a maximum of 5; at the published angles it is 0.185 with a maximum of 199. So the difference as shipped is the shader taking the sine and cosine of its own view angles at float32, which turns the camera by a small fraction of a degree and moves the silhouette by a pixel; it is not in the march. Where the two genuinely differ is sampling: the JS path shades a 200-cell grid and scales it up, the shader runs at the canvas backing resolution, and comparing them at their own native samplings mostly measures the upscale rather than the arithmetic, which is why the matched grid is the one quoted.

What the shader costs, on the record

A sphere trace at backing resolution is among the most expensive shaders in this atlas, and on a software renderer that is worth stating rather than discovering. Measured in the browser gate own Chrome, which on this machine has no GPU and runs ANGLE over SwiftShader, at a 318 by 256 card frame: this plate is 25.5 ms at the published constants and 30.2 at its most expensive power, against 42.6 for mandelbox, 21.9 for metamer and about 4 for mandelbrot, all four measured on the same harness in the same four consecutive runs and reproducible to a tenth of a millisecond. So it lands level with metamer, which held the record before mandelbox, and at three fifths of the plate it follows. The wave-21 pilot predicted this would be the worst frame in the atlas by a wide margin and it is not, and the reason is worth carrying: one estimate here costs 175.6 ns against 127 for mandelbox, but the object is compact and the estimate converges quickly on it, so the march is shorter — 8.9 estimates a ray at the card grid against 12.1 — and the two effects very nearly cancel. Four dials were traded and each was measured: the framing (2.9 of the containment radius rather than 2.6, 25.9 ms against 33.0), the obscurance (four directions by six steps rather than six by ten, mandelbox own measured trade), the shadow cone at 16 steps, and the march budget, which is not a lever here either. What cannot be traded is the resolution: a fragment shader cannot march at a capped grid the way the JS path does, because every fragment pays for its own ray. The reader who wants this plate cheap has GPU OFF, which puts it on the JS path.

Standing

Public domain. A published construction, posted openly in 2007 and 2009 and implemented here from the posted formulae rather than from the renderer of anybody at all — no source code from any fractal program was consulted or adapted. No patent has ever applied to it. The rendering method is published literature of the same standing: sphere tracing, the distance-bounded ray march for deterministic fractals, Lambert, Blinn and the ambient obscurance model are all papers, and none of them is or ever was encumbered.

Source

https://ff1.fractalforums.org/index.php?topic=742.0

TOUCHDESIGNER · GLSL

The same shader this plate runs, reframed for a GLSL TOP. Pasted bare it renders the published constants as a still frame; wire absTime.seconds into u_t on the Vectors page to animate it.

// FORMA — PL. 148 · MANDELBULB — Daniel White & Paul Nylander, 2007–2009 · after Rudy Rucker, 1988
//   vₖ₊₁ = vₖⁿ + c,   v₀ = 0,   c ∈ ℝ³
//   vⁿ = rⁿ(cos nθ·cos nφ,  sin nθ·cos nφ,  −sin nφ),   r = |v|
//   θ = atan2(y, x),   φ = atan(z / √(x²+y²)),   n = 8
//   dr → n·rⁿ⁻¹·dr + 1,   d = r·ln r / (2·dr)   — heuristic, not a bound
//   sphere trace: t → t + 0.75·d(o + t·û),   n̂ = ∇d / |∇d|
// TouchDesigner port — paste into a GLSL TOP's pixel shader. Set the
// resolution on the TOP's Common page. As pasted it renders the published
// constants as a still frame; to animate, add a uniform named u_t on the
// GLSL TOP's Vectors 1 page with the expression absTime.seconds.
// Constants are consts — edit to tweak; comments give the measured range.
// Written from the published mathematics, not adapted from any code.

#define u_res (uTDOutputInfo.res.zw)
uniform float u_t;               // absTime.seconds on the Vectors page; unset = still

const float u_phase = 0.3758;    // this plate's own grid phase, 0..1
// FORMA's FRACTALS accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.1389, 0.1912);
const vec3 u_pal_b = vec3(0.5, 0.151, 0.2078);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_power   = 8.0;         // n — power of the triplex map · live 2 .. 12
const float p_elev    = 22.0;        // view elevation ° · live -75 .. 75
const float p_az      = 28.0;        // view azimuth ° · live -180 .. 180
const float p_lightaz = -36.0;       // key light azimuth ° · live -180 .. 180
const float p_iter    = 12.0;        // iteration ceiling · live 4 .. 16
const float p_steps   = 160.0;       // march step budget · live 48 .. 192

/* The order's ramp — the same cosine formulation the JS kit uses, so a
   plate keeps its classification colour in either language. */
vec3 ramp(float t){
  return clamp(u_pal_a + u_pal_b * cos(6.28318530718 * (u_pal_c * t + u_pal_d)), 0.0, 1.0);
}

/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return fract(t / period + u_phase); }
float pingpong(float t, float period){
  float u = cycle(t, period);
  return u < 0.5 ? u * 2.0 : 2.0 - u * 2.0;
}


// The same march, at the canvas own backing resolution and in one pass.
// Every helper is bulb_ prefixed and nothing is declared at file scope, so
// a composition of this plate with itself still links.
float bulb_de(vec3 pos, out float trap){
  float n = floor(p_power + 0.5);
  float it = floor(p_iter + 0.5);
  vec3 v = pos;
  float dr = 1.0;
  float r = length(v);
  float mn = r;
  // 16 is the declared ceiling of p_iter -- mandelbrot fashion: a constant
  // loop bound at the slider own maximum, breaking on the live value.
  for (int i = 0; i < 16; i++){
    if (float(i) >= it) break;
    if (r > 4.0) break;
    float rho = sqrt(v.x * v.x + v.y * v.y);
    // The two degenerate directions, handled without an epsilon: on the
    // vertical axis the azimuth is 0 by convention, and at r = 0 the radius
    // factor is 0 so the direction cannot matter.
    float axis = rho > 0.0 ? 0.0 : 1.0;
    float orig = r > 0.0 ? 0.0 : 1.0;
    float azx = (v.x + axis) / (rho + axis);
    float azy = v.y / (rho + axis);
    float elx = rho / (r + orig);
    float ely = v.z / (r + orig);
    // de Moivre: n complex multiplies give cos(n.theta), sin(n.theta) and
    // cos(n.phi), sin(n.phi) with no trigonometry at all, and r^n alongside.
    // 12 is p_power's own ceiling.
    float cn = 1.0;
    float sn = 0.0;
    float cm = 1.0;
    float sm = 0.0;
    float rn = 1.0;
    for (int j = 0; j < 12; j++){
      if (float(j) >= n) break;
      float pa = cn * azx - sn * azy;
      sn = cn * azy + sn * azx;
      cn = pa;
      float pb = cm * elx - sm * ely;
      sm = cm * ely + sm * elx;
      cm = pb;
      rn *= r;
    }
    dr = n * (rn / (r + orig)) * dr + 1.0;
    v = vec3(rn * cn * cm, rn * sn * cm, -rn * sm) + pos;
    r = length(v);
    mn = min(mn, r);
  }
  trap = mn;
  return 0.5 * r * log(max(r, 1.0000001)) / dr;
}
float bulb_bound(float n){ return pow(2.0, 1.0 / (n - 1.0)); }
float bulb_eps(float tt, float pxa){ return max(2.0e-6, 0.30 * tt * pxa); }
float bulb_march(vec3 ro, vec3 rd, float t0, float t1, float pxa, out float trap){
  float tt = t0, hit = -1.0;
  trap = 0.0;
  for (int i = 0; i < 192; i++){
    if (float(i) >= p_steps) break;
    float tr;
    float d = bulb_de(ro + rd * tt, tr);
    if (d < bulb_eps(tt, pxa)){ trap = tr; hit = tt; break; }
    tt += d * 0.75;
    if (tt > t1) break;
  }
  return hit;
}
vec3 bulb_normal(vec3 pos, float e){
  float tr;
  float gx = bulb_de(pos + vec3(e, 0.0, 0.0), tr) - bulb_de(pos - vec3(e, 0.0, 0.0), tr);
  float gy = bulb_de(pos + vec3(0.0, e, 0.0), tr) - bulb_de(pos - vec3(0.0, e, 0.0), tr);
  float gz = bulb_de(pos + vec3(0.0, 0.0, e), tr) - bulb_de(pos - vec3(0.0, 0.0, e), tr);
  return normalize(vec3(gx, gy, gz));
}
float bulb_shadow(vec3 pos, vec3 l, float t0, float t1){
  float k = 1.0, tt = t0;
  for (int i = 0; i < 16; i++){
    float tr;
    float d = bulb_de(pos + l * tt, tr);
    if (d < 2.0e-6){ k = 0.0; break; }
    k = min(k, 11.0 * d / tt);
    tt += d * 0.75;
    if (tt > t1) break;
  }
  return clamp(k, 0.0, 1.0);
}
float bulb_obscurance(vec3 pos, vec3 n, float e, float R){
  vec3 u0 = vec3(0.0, 0.0, 1.0);
  if (abs(n.z) > 0.9) u0 = vec3(1.0, 0.0, 0.0);
  vec3 tg = normalize(cross(n, u0));
  vec3 bt = cross(n, tg);
  float sum = 0.0;
  for (int i = 0; i < 4; i++){
    float fi = (float(i) + 0.5) / 4.0;
    float rr = sqrt(fi), ph = float(i) * 2.39996323;
    vec3 d = normalize(tg * (rr * cos(ph)) + bt * (rr * sin(ph)) + n * sqrt(max(0.0, 1.0 - fi)));
    float tt = 2.0 * e, hit = R;
    for (int q = 0; q < 6; q++){
      float tr;
      float dd = bulb_de(pos + d * tt, tr);
      if (dd < e){ hit = tt; break; }
      tt += max(dd, e);
      if (tt >= R){ hit = R; break; }
    }
    sum += min(1.0, hit / R);
  }
  return sum / 4.0;
}
vec3 plate(vec2 uv){
  float DEG = 0.017453292519943295;
  float htan = 0.4142135623730951;   // tan of half a 45 degree field
  float n = floor(p_power + 0.5);
  float bound = bulb_bound(n);
  float dist = bound * 2.9;
  float az = p_az * DEG, el = p_elev * DEG;
  vec3 ro = vec3(dist * cos(el) * cos(az), dist * cos(el) * sin(az), dist * sin(el));
  vec3 fw = -ro / dist;
  // The world up axis is the map own polar axis, so the poles the theory
  // warns about are the top and bottom of the frame.
  vec3 rt = normalize(vec3(fw.y, -fw.x, 0.0));
  vec3 up = cross(rt, fw);

  float la = p_lightaz * DEG;
  vec3 kc = vec3(sin(la) * 0.7762, 0.6305, -cos(la) * 0.7762);
  vec3 key = rt * kc.x + up * kc.y + fw * kc.z;
  vec3 fil = rt * 0.9048 + up * -0.2088 + fw * -0.3712;

  float asp = u_res.x / u_res.y;
  float pxa = 2.0 * htan / u_res.y;
  float sx = uv.x * 2.0 - 1.0, sy = 1.0 - uv.y * 2.0;
  vec3 rd = normalize(fw + rt * (sx * htan * asp) + up * (sy * htan));

  vec3 col = vec3(4.0, 6.0, 10.0) / 255.0;
  float bq = dot(ro, rd);
  float cq = dot(ro, ro) - bound * bound;
  float disc = bq * bq - cq;
  if (disc > 0.0){
    float sq = sqrt(disc);
    float t0 = max(0.0, -bq - sq), t1 = -bq + sq;
    if (t1 > 0.0){
      float trap;
      float th = bulb_march(ro, rd, t0, t1, pxa, trap);
      if (th > 0.0){
        vec3 pos = ro + rd * th;
        float e = bulb_eps(th, pxa);
        vec3 nm = bulb_normal(pos, 1.2 * e);
        float ndl = max(0.0, dot(nm, key));
        float sh = 0.0;
        if (ndl > 0.0) sh = bulb_shadow(pos + nm * (3.0 * e), key, 3.0 * e, bound * 2.0);
        float ndf = max(0.0, dot(nm, fil));
        float ao = bulb_obscurance(pos, nm, e, 0.09 * bound);
        vec3 hv = normalize(key - rd);
        float spec = pow(max(0.0, dot(nm, hv)), 36.0) * ndl * sh;
        // trap^(n-1)/2 is the closest approach measured on the escape
        // criterion own scale, so the colouring means the same thing at
        // every power rather than needing a window measured at one.
        vec3 al = ramp(0.82 + 0.30 * clamp(pow(trap, n - 1.0) * 0.5, 0.0, 1.0));
        col = al * (0.86 * ndl * sh + 0.18 * ndf)
            + al * vec3(0.4, 0.6, 1.0) * (0.25 * ao)
            + vec3(1.0) * (0.30 * spec);
      }
    }
  }
  return col;
}

out vec4 fragColor;
void main(){
  // FORMA's uv runs y-down, matching its canvas; TD's vUV runs up
  vec2 uv = vec2(vUV.s, 1.0 - vUV.t);
  fragColor = TDOutputSwizzle(vec4(plate(uv), 1.0));
}

NUKE · BLINKSCRIPT

The same shader this plate runs, transpiled to a BlinkScript kernel. Paste it into a BlinkScript node's Kernel Source and press Recompile; every constant arrives as a knob at its published value, and u_t animates with the expression frame/24. Compiled and rendered in Nuke 17.1, then compared against this plate on the page.

// FORMA — PL. 148 · MANDELBULB — Daniel White & Paul Nylander, 2007–2009 · after Rudy Rucker, 1988
//   vₖ₊₁ = vₖⁿ + c,   v₀ = 0,   c ∈ ℝ³
//   vⁿ = rⁿ(cos nθ·cos nφ,  sin nθ·cos nφ,  −sin nφ),   r = |v|
//   θ = atan2(y, x),   φ = atan(z / √(x²+y²)),   n = 8
//   dr → n·rⁿ⁻¹·dr + 1,   d = r·ln r / (2·dr)   — heuristic, not a bound
//   sphere trace: t → t + 0.75·d(o + t·û),   n̂ = ∇d / |∇d|
// Nuke port — a BlinkScript kernel. Paste into a BlinkScript node's Kernel
// Source and press Recompile. Every constant arrives as a knob at its published
// value (the comment gives the measured range); u_t is a knob too — animate it
// with the expression frame/24 or leave it at 0 for the still frame. Written
// from the published mathematics, not adapted from any code.
// Transpiled from the shader this plate runs on the page (GLSL ES 3.00):
// vec → float2/3/4, swizzles expanded, GLSL builtins Blink lacks written out
// as forma_ functions, float literals suffixed. Compiled and rendered in a
// real Nuke (17.1v1) and compared against this plate on the page: 34 of 34.
//
// plate() and its helpers are written to a single exit — the loop that runs
// once. That is not a style: Blink 17.1 drops a conditional early return from
// a called function while Vectorize is on, which is the node default, with no
// warning and no error. Written this way it paints correctly as pasted.

kernel Forma_mandelbulb : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_power;   // n — power of the triplex map · live 2 .. 12
  float p_elev;    // view elevation ° · live -75 .. 75
  float p_az;      // view azimuth ° · live -180 .. 180
  float p_lightaz; // key light azimuth ° · live -180 .. 180
  float p_iter;    // iteration ceiling · live 4 .. 16
  float p_steps;   // march step budget · live 48 .. 192

local:
  float2 u_res;
  float u_phase;
  float3 u_pal_a, u_pal_b, u_pal_c, u_pal_d;

  void define(){
    defineParam(u_t, "u_t", 0.0f);
    defineParam(p_power, "p_power", 8.0f);
    defineParam(p_elev, "p_elev", 22.0f);
    defineParam(p_az, "p_az", 28.0f);
    defineParam(p_lightaz, "p_lightaz", -36.0f);
    defineParam(p_iter, "p_iter", 12.0f);
    defineParam(p_steps, "p_steps", 160.0f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.3758f;    // this plate's own grid phase, 0..1
    // FORMA's FRACTALS accent as cosine-gradient coefficients
    u_pal_a = float3(0.46f, 0.1389f, 0.1912f);
    u_pal_b = float3(0.5f, 0.151f, 0.2078f);
    u_pal_c = float3(1.0f, 1.0f, 1.0f);
    u_pal_d = float3(0.0f, 0.05f, 0.1f);
  }

  /* GLSL builtins Blink lacks, written as templates rather than overload sets.
     Blink's operators return expression templates (Swizzle<float,N>), so a call
     passing an expression cannot resolve against an overload set on float2
     against float3 — measured in Nuke 17.1: a float2 expression is ambiguous
     between the two, while scalar-against-vector resolves. A template deduces
     the expression's own type, so the ambiguity cannot arise. */
  template <class T> T forma_fract(T v){ return v - floor(v); }
  template <class T, class S> T forma_mod(T x, S y){ return x - y * floor(x / y); }
  /* Blink's own min/max/clamp take no scalar bound against a vector, which GLSL
     does; v * 0.0f + b is that bound at the vector's own width, and collapses to
     b when v is a scalar, so one template serves both. */
  template <class T, class S> T forma_min(T a, S b){ return min(a, a * 0.0f + b); }
  template <class T, class S> T forma_max(T a, S b){ return max(a, a * 0.0f + b); }
  template <class T, class S> T forma_clamp(T v, S lo, S hi){ return clamp(v, v * 0.0f + lo, v * 0.0f + hi); }
  int forma_min(int a, int b){ return min(a, b); }
  int forma_max(int a, int b){ return max(a, b); }
  /* GLSL step(edge, x) is 1 where x >= edge; floor(sign(x - e) * 0.5 + 1) is
     that exactly, equality included, out of builtins Blink does have. */
  template <class T, class S> T forma_step(S e, T x){ return floor(sign(x - e) * 0.5f + 1.0f); }
  template <class T, class S> T forma_smoothstep(S a, S b, T x){
    T t = forma_clamp((x - a) / (b - a), 0.0f, 1.0f);
    return t * t * (3.0f - 2.0f * t);
  }
  template <class T> float forma_distance(T a, T b){ return length(a - b); }
  float forma_tanh(float x){ float e = exp(2.0f * x); return (e - 1.0f) / (e + 1.0f); }
  float forma_radians(float d){ return d * 0.01745329252f; }
  // the page's hash2 is exact uint32; Blink has int, so the shifts are made
  // logical by masking and the read-back is lifted into 0 .. 2^32
  /* A uint read back as a float. Blink has no unsigned type, so a value past
     2^31 arrives as a negative int and float() of it is negative. Measured on
     gabor, whose own generator then returned uniforms in [-0.5, 0.5) and drew
     a different picture — it compiled, it rendered, and only comparing it with

  /* The order's ramp — the same cosine formulation the JS kit uses, so a
     plate keeps its classification colour in either language. */
  float3 ramp(float t){
    return forma_clamp(u_pal_a + u_pal_b * cos(6.28318530718f * (u_pal_c * t + u_pal_d)), 0.0f, 1.0f);
  }

  /* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
  float cycle(float t, float period){ return forma_fract(t / period + u_phase); }
  float pingpong(float t, float period){
    float u = cycle(t, period);
    return u < 0.5f ? u * 2.0f : 2.0f - u * 2.0f;
  }


  // The same march, at the canvas own backing resolution and in one pass.
  // Every helper is bulb_ prefixed and nothing is declared at file scope, so
  // a composition of this plate with itself still links.
  float bulb_de(float3 pos, float& trap){
    float n = floor(p_power + 0.5f);
    float it = floor(p_iter + 0.5f);
    float3 v = pos;
    float dr = 1.0f;
    float r = length(v);
    float mn = r;
    // 16 is the declared ceiling of p_iter -- mandelbrot fashion: a constant
    // loop bound at the slider own maximum, breaking on the live value.
    for (int i = 0; i < 16; i++){
      if (float(i) >= it) break;
      if (r > 4.0f) break;
      float rho = sqrt(v.x * v.x + v.y * v.y);
      // The two degenerate directions, handled without an epsilon: on the
      // vertical axis the azimuth is 0 by convention, and at r = 0 the radius
      // factor is 0 so the direction cannot matter.
      float axis = rho > 0.0f ? 0.0f : 1.0f;
      float orig = r > 0.0f ? 0.0f : 1.0f;
      float azx = (v.x + axis) / (rho + axis);
      float azy = v.y / (rho + axis);
      float elx = rho / (r + orig);
      float ely = v.z / (r + orig);
      // de Moivre: n complex multiplies give cos(n.theta), sin(n.theta) and
      // cos(n.phi), sin(n.phi) with no trigonometry at all, and r^n alongside.
      // 12 is p_power's own ceiling.
      float cn = 1.0f;
      float sn = 0.0f;
      float cm = 1.0f;
      float sm = 0.0f;
      float rn = 1.0f;
      for (int j = 0; j < 12; j++){
        if (float(j) >= n) break;
        float pa = cn * azx - sn * azy;
        sn = cn * azy + sn * azx;
        cn = pa;
        float pb = cm * elx - sm * ely;
        sm = cm * ely + sm * elx;
        cm = pb;
        rn *= r;
      }
      dr = n * (rn / (r + orig)) * dr + 1.0f;
      v = float3(rn * cn * cm, rn * sn * cm, -rn * sm) + pos;
      r = length(v);
      mn = forma_min(mn, r);
    }
    trap = mn;
    return 0.5f * r * log(forma_max(r, 1.0000001f)) / dr;
  }
  float bulb_bound(float n){ return pow(2.0f, 1.0f / (n - 1.0f)); }
  float bulb_eps(float tt, float pxa){ return forma_max(2.0e-6f, 0.30f * tt * pxa); }
  float bulb_march(float3 ro, float3 rd, float t0, float t1, float pxa, float& trap){
    float tt = t0;
    float hit = -1.0f;
    trap = 0.0f;
    for (int i = 0; i < 192; i++){
      if (float(i) >= p_steps) break;
      float tr;
      float d = bulb_de(ro + rd * tt, tr);
      if (d < bulb_eps(tt, pxa)){ trap = tr; hit = tt; break; }
      tt += d * 0.75f;
      if (tt > t1) break;
    }
    return hit;
  }
  float3 bulb_normal(float3 pos, float e){
    float tr;
    float gx = bulb_de(pos + float3(e, 0.0f, 0.0f), tr) - bulb_de(pos - float3(e, 0.0f, 0.0f), tr);
    float gy = bulb_de(pos + float3(0.0f, e, 0.0f), tr) - bulb_de(pos - float3(0.0f, e, 0.0f), tr);
    float gz = bulb_de(pos + float3(0.0f, 0.0f, e), tr) - bulb_de(pos - float3(0.0f, 0.0f, e), tr);
    return normalize(float3(gx, gy, gz));
  }
  float bulb_shadow(float3 pos, float3 l, float t0, float t1){
    float k = 1.0f;
    float tt = t0;
    for (int i = 0; i < 16; i++){
      float tr;
      float d = bulb_de(pos + l * tt, tr);
      if (d < 2.0e-6f){ k = 0.0f; break; }
      k = forma_min(k, 11.0f * d / tt);
      tt += d * 0.75f;
      if (tt > t1) break;
    }
    return forma_clamp(k, 0.0f, 1.0f);
  }
  float bulb_obscurance(float3 pos, float3 n, float e, float R){
    float3 u0 = float3(0.0f, 0.0f, 1.0f);
    if (fabs(n.z) > 0.9f) u0 = float3(1.0f, 0.0f, 0.0f);
    float3 tg = normalize(cross(n, u0));
    float3 bt = cross(n, tg);
    float sum = 0.0f;
    for (int i = 0; i < 4; i++){
      float fi = (float(i) + 0.5f) / 4.0f;
      float rr = sqrt(fi);
      float ph = float(i) * 2.39996323f;
      float3 d = normalize(tg * (rr * cos(ph)) + bt * (rr * sin(ph)) + n * sqrt(forma_max(0.0f, 1.0f - fi)));
      float tt = 2.0f * e;
      float hit = R;
      for (int q = 0; q < 6; q++){
        float tr;
        float dd = bulb_de(pos + d * tt, tr);
        if (dd < e){ hit = tt; break; }
        tt += forma_max(dd, e);
        if (tt >= R){ hit = R; break; }
      }
      sum += forma_min(1.0f, hit / R);
    }
    return sum / 4.0f;
  }
  float3 plate(float2 uv){
    float DEG = 0.017453292519943295f;
    float htan = 0.4142135623730951f;   // tan of half a 45 degree field
    float n = floor(p_power + 0.5f);
    float bound = bulb_bound(n);
    float dist = bound * 2.9f;
    float az = p_az * DEG;
    float el = p_elev * DEG;
    float3 ro = float3(dist * cos(el) * cos(az), dist * cos(el) * sin(az), dist * sin(el));
    float3 fw = -ro / dist;
    // The world up axis is the map own polar axis, so the poles the theory
    // warns about are the top and bottom of the frame.
    float3 rt = normalize(float3(fw.y, -fw.x, 0.0f));
    float3 up = cross(rt, fw);

    float la = p_lightaz * DEG;
    float3 kc = float3(sin(la) * 0.7762f, 0.6305f, -cos(la) * 0.7762f);
    float3 key = rt * kc.x + up * kc.y + fw * kc.z;
    float3 fil = rt * 0.9048f + up * -0.2088f + fw * -0.3712f;

    float asp = u_res.x / u_res.y;
    float pxa = 2.0f * htan / u_res.y;
    float sx = uv.x * 2.0f - 1.0f;
    float sy = 1.0f - uv.y * 2.0f;
    float3 rd = normalize(fw + rt * (sx * htan * asp) + up * (sy * htan));

    float3 col = float3(4.0f, 6.0f, 10.0f) / 255.0f;
    float bq = dot(ro, rd);
    float cq = dot(ro, ro) - bound * bound;
    float disc = bq * bq - cq;
    if (disc > 0.0f){
      float sq = sqrt(disc);
      float t0 = forma_max(0.0f, -bq - sq);
      float t1 = -bq + sq;
      if (t1 > 0.0f){
        float trap;
        float th = bulb_march(ro, rd, t0, t1, pxa, trap);
        if (th > 0.0f){
          float3 pos = ro + rd * th;
          float e = bulb_eps(th, pxa);
          float3 nm = bulb_normal(pos, 1.2f * e);
          float ndl = forma_max(0.0f, dot(nm, key));
          float sh = 0.0f;
          if (ndl > 0.0f) sh = bulb_shadow(pos + nm * (3.0f * e), key, 3.0f * e, bound * 2.0f);
          float ndf = forma_max(0.0f, dot(nm, fil));
          float ao = bulb_obscurance(pos, nm, e, 0.09f * bound);
          float3 hv = normalize(key - rd);
          float spec = pow(forma_max(0.0f, dot(nm, hv)), 36.0f) * ndl * sh;
          // trap^(n-1)/2 is the closest approach measured on the escape
          // criterion own scale, so the colouring means the same thing at
          // every power rather than needing a window measured at one.
          float3 al = ramp(0.82f + 0.30f * forma_clamp(pow(trap, n - 1.0f) * 0.5f, 0.0f, 1.0f));
          col = al * (0.86f * ndl * sh + 0.18f * ndf)
              + al * float3(0.4f, 0.6f, 1.0f) * (0.25f * ao)
              + float3(1.0f) * (0.30f * spec);
        }
      }
    }
    return col;
  }

  void process(int2 pos){
    // FORMA's uv runs y-down like its canvas; Nuke's rows run up
    float2 uv = float2((float(pos.x) + 0.5f) / u_res.x, 1.0f - (float(pos.y) + 0.5f) / u_res.y);
    float3 c = plate(uv);
    dst() = float4(c.x, c.y, c.z, 1.0f);
  }
};