Mixwell.hlsl
The complete Mixwell library: brushes, particle drift, Reverse-Drift Functions, and the Newton solver.
// SPDX-License-Identifier: MIT
// Copyright (c) 2026 Doug L. James and Ethan James
///////////////////////////////////////////////////////////////////////
//// Mixwell Library Functions /// @author Doug L James, 2026 ////////
///////////////////////////////////////////////////////////////////////
static const float PI = 3.14159265358979323846f; // π
static const float PI_2 = 1.57079632679489661923f; // π/2
// Cutoff distance (as multiple of tine radius) beyond which rdSegment is zero.
static const float RDSEGMENT_MASK_R_FACTOR = 10.0f;
// Euclidean mod // NEW
float emod(float x, float n)
{
return frac(x / n) * n;
}
float2 emod(float2 x, float2 n)
{
return frac(x / n) * n;
}
float3 emod(float3 x, float3 n)
{
return frac(x / n) * n;
}
/*** [BEGIN] SDF FUNCTIONS *******************************************/
float sdCircle(float2 p, float r)
{
return length(p) - r;
}
float sdBox(float2 p, float2 b)
{
float2 d = abs(p) - b;
return length(max(d, float2(0.0f, 0.0f))) + min(max(d.x, d.y), 0.0f);
}
// Unsigned distance from the infinite line specified by a ray origin (ro) and direction (unit rd).
float udLine(float2 p, float2 ro, float2 rd)
{
float2 v = p - ro;
return length(v - dot(v, rd) * rd / dot(rd, rd));
}
float sdSegment(float2 p, float2 a, float2 b)
{
float2 pa = p - a;
float2 ba = b - a;
float h = clamp(dot(pa, ba) / dot(ba, ba), 0.0f, 1.0f);
return length(pa - ba * h);
}
float sdTriangle(float2 p, float2 p0, float2 p1, float2 p2)
{
float2 e0 = p1 - p0;
float2 e1 = p2 - p1;
float2 e2 = p0 - p2;
float2 v0 = p - p0;
float2 v1 = p - p1;
float2 v2 = p - p2;
float2 pq0 = v0 - e0 * clamp(dot(v0, e0) / dot(e0, e0), 0.0f, 1.0f);
float2 pq1 = v1 - e1 * clamp(dot(v1, e1) / dot(e1, e1), 0.0f, 1.0f);
float2 pq2 = v2 - e2 * clamp(dot(v2, e2) / dot(e2, e2), 0.0f, 1.0f);
float s = sign(e0.x * e2.y - e0.y * e2.x);
float2 d = min(min(float2(dot(pq0, pq0), s * (v0.x * e0.y - v0.y * e0.x)),
float2(dot(pq1, pq1), s * (v1.x * e1.y - v1.y * e1.x))),
float2(dot(pq2, pq2), s * (v2.x * e2.y - v2.y * e2.x)));
return -sqrt(d.x) * sign(d.y);
}
/*** [END] SDF FUNCTIONS *******************************************/
/*****************************************************************************/
/*** [BEGIN] COLOR & TEST IMAGES ***********************/
/*****************************************************************************/
// Matlab Jet colormap :/
float4 jet(float x)
{
float r = clamp((x < 0.7f) ? (4.0f * x - 1.5f) : (-4.0f * x + 4.5f), 0.0f, 1.0f);
float g = clamp((x < 0.5f) ? (4.0f * x - 0.5f) : (-4.0f * x + 3.5f), 0.0f, 1.0f);
float b = clamp((x < 0.3f) ? (4.0f * x + 0.5f) : (-4.0f * x + 2.5f), 0.0f, 1.0f);
return float4(r, g, b, 1.0f);
}
float3 rgb(int r, int g, int b)
{
return float3(float(r), float(g), float(b)) / 255.0f;
}
float3 quantizeColor(float3 c, float n, int bands)
{
float3 res = floor(c * n) / n;
if (bands)
res *= smoothstep(-0.0001f, +0.0001f, length(c - res) - 0.005f);
return res;
}
// Paint texture
float3 paintChocolate(float2 q)
{
float theta = 0.3f;
float c = cos(theta);
float s = sin(theta);
float2 p = float2(c * q.x - s * q.y, s * q.x + c * q.y);
float k0 = 45.0f;
float k1 = 29.0f;
float phi0 = 3.0f;
float phi1 = 2.0f;
float3 col0 = float3(0.6862745098f, 0.4f, 0.1490196078f); // rgb(175, 102, 38);
float3 col1 = float3(0.6901960784f, 0.9411764706f, 0.9921568627f); // rgb(176, 240, 253);
float3 colChocolate = float3(0.1490196078f, 0.0862745098f, 0.0901960784f); // rgb(38,22,23);
float splat0 = smoothstep(0.5f, 0.64f,
sin(k0 * p.x * (1.0f - 0.153f * p.x) + phi0) * sin(k0 * p.y * (1.0f + 0.6f * p.x) + phi0) *
(1.0f - 0.2f * sin(12.0f * p.x * p.x)));
float splat1 = smoothstep(0.5f, 0.64f,
sin(k1 * p.x * (1.0f + 0.230f * p.x) + phi1) * cos(k1 * p.y * (1.0f - 0.1f * p.x) + phi1) *
(1.0f - 0.2f * sin(12.0f * p.x * p.y)));
float3 col = lerp(colChocolate + splat0 * col0, col1, splat1);
return col;
}
float3 paintTestSquare(float2 q) // NEW
{
return float3(1.0f, 1.0f, 1.0f) * smoothstep(0.1f, 0.11f, q.x) * smoothstep(0.2f, 0.19f, q.x) *
smoothstep(0.1f, 0.11f, q.y) * smoothstep(0.2f, 0.19f, q.y);
}
float3 paintStripes(float2 q)
{
return pow(smoothstep(0.0f, 1.0f, frac(6.0f * q.y + sin(33.0f * q.x))), 3.0f) * float3(1.0f, 1.0f, 1.0f);
}
// Adapted from https://webgl-operate.org/examples/canvassize-example.html // NEW
float3 paintTestImage(float2 fragCoord)
{
const float CELL_WIDTH = 1.0f / 32.0f;
float3 x3 = float3(fragCoord.x, fragCoord.x + 1.0f, fragCoord.x + 2.0f);
float3 y3 = float3(fragCoord.y, fragCoord.y + 1.0f, fragCoord.y + 2.0f);
float3 x = step(emod(x3, 3.0f), 1.0f); // emod
float3 y = step(emod(y3, 3.0f), 1.0f); // emod
float cell = step(emod(fragCoord.x * CELL_WIDTH + floor(fragCoord.y * CELL_WIDTH), 2.0f), 1.0f); // emod
return lerp(x, y, cell);
}
float3 paintCheckerboard(float2 q, float h) // NEW
{
int2 cell = int2(floor(q / h));
return (cell.x + cell.y) & 1;
}
float3 splatBlobRows(float2 p, float3 colBG, float3 colBlob, float HX, float HY, float dY, float blobRadius, float isOddRow) // NEW
{
float Yoffset = (isOddRow) ? 0.5f * HY : 0.f;
// SNAP p to closest horizontal row:
float ySnap = Yoffset + round((p.y - Yoffset) / HY) * HY; // y of closest row
p.y -= ySnap;
p.x -= sin(437.5453 * ySnap) * 437.5453; // SHIFT X per ROW
// STRING OF BLOBS ON y=0 AXIS:
float F = 1239.f + ySnap; // nat freq for variation
// Mission: Find our blob i:
float x0 = HX * floor(p.x / HX);
float x1 = x0 + HX; // ceil
float2 p0 = float2(x0, dY * sin(F * x0));
float2 p1 = float2(x1, dY * sin(F * x1));
float D0 = dot(p - p0, p - p0);
float D1 = dot(p - p1, p - p1);
float iC = (D0 < D1) ? round(x0 / HX) : round(x1 / HX);
// Render blob: It takes three blobs to raise a blob:
float xC = iC * HX; // center (i)
float xL = xC - HX; // left (i-1)
float xR = xC + HX; // right (i+1)
float yC = dY * sin(F * xC); // ignore Yoffset here
float yL = dY * sin(F * xL);
float yR = dY * sin(F * xR);
float2 pC = float2(xC, yC);
float2 pL = float2(xL, yL);
float2 pR = float2(xR, yR);
float e2 = 0.00001 * HX * HX;
float DC = dot(p - pC, p - pC) + e2;
float DL = dot(p - pL, p - pL) + e2;
float DR = dot(p - pR, p - pR) + e2;
float f = 1.f / DC - 1.f / DL - 1.f / DR;
// Implicit blob size hint:
f -= 1.f / (blobRadius * blobRadius);
float reg = smoothstep(HX, 2.f * HX, max(0.f, f));
return lerp(colBG, colBlob, reg);
}
static const float3 paletteBluePurple[6] =// NEW
{
float3(8, 8, 140) / 255.0f,
float3(31, 17, 112) / 255.0f,
float3(17, 151, 247) / 255.0f,
float3(97, 25, 191) / 255.0f,
float3(235, 251, 252) / 255.0f,
float3(255, 212, 235) / 255.0f,
};
//------------------------------------------------------------------------------
// Palette 1: Coolors example (brown blue)
// 1a344d,1e8eb6,fc9843,5a1800,f7f5f6
//------------------------------------------------------------------------------
static const float3 paletteCoolors1[6] =// NEW
{
float3(0.1019608, 0.2039216, 0.3019608),
float3(0.3529412, 0.09411765, 0.0),
float3(0.1176471, 0.5568628, 0.7137255),
float3(0.9882353, 0.5960785, 0.2627451),
float3(0.9686275, 0.9607843, 0.9647059),
float3(0.9686275, 0.9607843, 0.9647059),
};
//------------------------------------------------------------------------------
// Palette: fb6107,f3de2c,7cb518,5c8001,fbb02d
// Names: Blaze Orange, Golden Glow, Lime Moss, Forest Moss, Sunflower Gold
//------------------------------------------------------------------------------
static const float3 paletteGoldenMeadow[6] =// NEW
{
float3(0.9843137, 0.3803922, 0.0274510), // fb6107
float3(0.3607843, 0.5019608, 0.0039216), // 5c8001
float3(0.9843137, 0.6901961, 0.1764706), // fbb02d
float3(0.4862745, 0.7098039, 0.0941176), // 7cb518
float3(0.9529412, 0.8705882, 0.1725490), // f3de2c
float3(0.9529412, 0.8705882, 0.1725490), // f3de2c
};
float3 paintBlobsWithPalette(float2 p, const float3 palette[6])// NEW
{
float HX = 0.25; // X CELL SIZE
float HY = 1.f; // ROW SPACING
float dY = HX / 10.f; // Y POSITION VARIATION
float blobRadius = 2.f * HX; // BLOB SIZE HINT FOR SMALLER BLOBS
float3 col = float3(0, 0, 0); // Background
col = splatBlobRows(p, col, palette[0], HX, HY, dY, 2.f * HX, false);
col = splatBlobRows(p, col, palette[1], HX, HY, dY, 2.f * HX, true);
col = splatBlobRows(p, col, palette[2], HX, HY, dY, HX / 2., false);
col = splatBlobRows(p, col, palette[3], HX, HY, dY, HX / 3., true);
col = splatBlobRows(p, col, palette[4], HX, HY, dY, HX / 6., false);
col = splatBlobRows(p, col, palette[5], HX, HY, dY, HX / 7., true);
return col;
}
float3 paintBlobs(float2 p) // NEW
{
return paintBlobsWithPalette(p, paletteGoldenMeadow);
}
/*****************************************************************************/
/*** [ END ] COLOR & TEST IMAGES ***********************/
/*****************************************************************************/
/*** [BEGIN] UTILITIES *******************************************/
float2 rotate(float2 v, float theta)
{
float c = cos(theta);
float s = sin(theta);
return float2(c * v.x - s * v.y, s * v.x + c * v.y);
}
// Returns (-y,x) given v=(x,y).
float2 rot90(float2 v)
{
return float2(-v.y, v.x);
}
// Normalizes p if ‖p‖<1.
float2 projectOutsideUnitDisk(float2 p)
{
float R2 = p.x * p.x + p.y * p.y;
return (R2 >= 1.f) ? p : p / sqrt(R2);
}
// Simple pseudorandom float-to-float[0,1] hash.
float iqhash(float n)
{
return frac(sin(n) * 43758.5453);
}
/*** [END] UTILITIES *******************************************/
/*****************************************************************************/
/*** [BEGIN] CARLSON ITERATIONS AND ELLIPTIC FUNCTIONS ***********************/
/*****************************************************************************/
#define CARLSON_MAX_ITER 15 // Maximum number of iterations
#define CARLSON_TOL 0.001f // Tolerance for convergence
// Carlson's RF function [Carlson 1995, p.16] (as .x) and the final "n" iteration# (as .y).
float2 CarlsonRFn(float x, float y, float z)
{
float A0 = (x + y + z) * 0.3333333333333333f;
float An = A0;
float Q = pow(3.0f * CARLSON_TOL, -1.0f / 6.0f) * max(max(abs(A0 - x), abs(A0 - y)), abs(A0 - z));
float xn = x;
float yn = y;
float zn = z;
float n4 = 1.0f;
int n;
for (n = 0; n < CARLSON_MAX_ITER; n++)
{
if (Q <= CARLSON_TOL * An)
{
break;
}
float lambda = sqrt(xn * yn) + sqrt(xn * zn) + sqrt(yn * zn);
An = (An + lambda) * 0.25f;
xn = (xn + lambda) * 0.25f;
yn = (yn + lambda) * 0.25f;
zn = (zn + lambda) * 0.25f;
Q *= 0.25f;
n4 *= 4.0f;
}
float X = (A0 - x) / (n4 * An);
float Y = (A0 - y) / (n4 * An);
float Z = -X - Y;
float E2 = X * Y - Z * Z;
float E3 = X * Y * Z;
float RF = (1.f - (E2 / 10.f) + (E3 / 14.f) + (E2 * E2 / 24.f) - (3.f / 44.f * E2 * E3)) / sqrt(An);
return float2(RF, float(n));
}
// Carlson's RD function [Carlson 1995, p.20] (as .x) and the final "n" iteration# (as .y).
float2 CarlsonRDn(float x, float y, float z)
{
float A0 = (x + y + 3.f * z) / 5.f;
float An = A0;
float Q = pow(CARLSON_TOL / 4.f, -1.f / 6.f) * max(max(abs(A0 - x), abs(A0 - y)), abs(A0 - z));
float xn = x;
float yn = y;
float zn = z;
float Sum = 0.f;
float n4 = 1.f;
int n;
for (n = 0; n < CARLSON_MAX_ITER; n++)
{
float lambda = sqrt(xn * yn) + sqrt(xn * zn) + sqrt(yn * zn);
Sum += 1.f / (n4 * sqrt(zn) * (zn + lambda));
if (Q <= CARLSON_TOL * An)
{
break;
}
An = (An + lambda) * 0.25f;
xn = (xn + lambda) * 0.25f;
yn = (yn + lambda) * 0.25f;
zn = (zn + lambda) * 0.25f;
Q *= 0.25f;
n4 *= 4.f;
}
float X = (A0 - x) / (n4 * An);
float Y = (A0 - y) / (n4 * An);
float Z = -(X + Y) / 3.f;
float XY = X * Y;
float Z2 = Z * Z;
float E2 = XY - 6.f * Z2;
float E3 = (3.f * XY - 8.f * Z2) * Z;
float E4 = 3.f * (XY - Z2) * Z2;
float E5 = XY * Z2 * Z;
float Expansion = (1.f - 3.f / 14.f * E2 + E3 / 6.f + 9.f / 88.f * E2 * E2 - 3.f / 22.f * E4 - 9.f / 52.f * E2 * E3 + 3.f / 26.f * E5);
float RD = Expansion / (n4 * An * sqrt(An)) + 3.f * Sum;
return float2(RD, float(n));
}
// Carlson's RF function [Carlson 1995, p.16]
float CarlsonRF(float x, float y, float z)
{
return CarlsonRFn(x, y, z).x;
}
// Carlson's RD function [Carlson 1995, p.20]
float CarlsonRD(float x, float y, float z)
{
return CarlsonRDn(x, y, z).x;
}
// Mathematica-style Incomplete Elliptic Integral of the First Kind, EllipticF(φ,m), with m=k² parameter.
// Computed in terms of Carlson's RF function [Carlson 1995, (4.5)].
float EllipticF(float phi, float m)
{
float sinPhi = sin(phi);
float cosPhi = cos(phi);
return sinPhi * CarlsonRF(cosPhi * cosPhi, 1.f - m * sinPhi * sinPhi, 1.f); /// 1.0 vs cosPhi
}
// Mathematica-style Incomplete Elliptic Integral of the Second Kind, EllipticE(φ,m), with m=k² parameter.
// Computed in terms of Carlson's RF and RD functions [Carlson 1995, (4.6)].
float EllipticE(float phi, float m)
{
float sinPhi = sin(phi);
float cosPhi = cos(phi);
float cos2Phi = cosPhi * cosPhi;
float sin2Phi = sinPhi * sinPhi;
float sin3Phi = sin2Phi * sinPhi;
return sinPhi * CarlsonRF(cos2Phi, 1.f - m * sin2Phi, 1.f) -
(m * sin3Phi / 3.f) * CarlsonRD(cos2Phi, 1.f - m * sin2Phi, 1.f);
}
// Mathematica-style Incomplete Elliptic Integrals of both the First and Second Kinds,
// EllipticF(φ,m) and EllipticE(φ,m), respectively, with m=k² parameter.
// Faster than independent calls since shares a common CarlsonRF evaluation.
// Also returns the max number of Carlson RD | RF iterations, n.
// @return float3(F,E,n)
float3 EllipticFEn(float phi, float m)
{
float sinPhi = sin(phi);
float cosPhi = cos(phi);
float sinPhi2 = sinPhi * sinPhi;
float sinPhi3 = sinPhi2 * sinPhi;
float x = cosPhi * cosPhi;
float y = 1.f - m * sinPhi2;
float2 RFn = CarlsonRFn(x, y, 1.f); // (RF,nRF)
float2 RDn = CarlsonRDn(x, y, 1.f); // (RD,nRD)
float n = max(RFn.y, RDn.y); // max(nRF,nRD)
float F = sinPhi * RFn.x;
float E = F - (m * sinPhi3 / 3.f) * RDn.x;
return float3(F, E, n);
}
// Mathematica-style Incomplete Elliptic Integrals of both the First and Second Kinds,
// EllipticF(φ,m) and EllipticE(φ,m), respectively, with m=k² parameter.
// Faster than independent calls since shares a common CarlsonRF evaluation.
// @return float2(F,E)
float2 EllipticFE(float phi, float m)
{
return EllipticFEn(phi, m).xy;
}
/*****************************************************************************/
/*** [END] CARLSON ITERATIONS AND ELLIPTIC FUNCTIONS *************************/
/*****************************************************************************/
/*****************************************************************************/
/*** [BEGIN] MAXWELL FLOW PROBLEM AND MONOTONE-PRECONDITIONED NEWTON SOLVER **/
/*****************************************************************************/
// Clamps |σ| ≤ (1-ε)π/2.
float mxwClampSigmaPosY(float sigma)
{
float sigmaMax = 0.99999 * PI_2; // avoid singularity
return clamp(sigma, -sigmaMax, sigmaMax);
}
// Gradient of η_∞(p)
float2 mxwGradEtaInf(float2 p)
{
float X = p.x;
float Y = p.y;
float R2 = dot(p, p);
float R4 = R2 * R2;
return float2(2.f * X * Y / R4, 1.f + 2.f * Y * Y / R4 - 1.f / R2);
}
// Signed η_∞ for point p=(X,Y) outside the unit circle.
float mxwEtaInf(float2 p)
{
float X = p.x;
float Y = p.y;
return Y - Y / (X * X + Y * Y);
}
// Maxwell-flow angle σ((X,Y))∈(π,-π), where σ=0 corresponds to the +y axis.
// Input: p - Point outside the unit circle.
float mxwSigmaAt(float2 p)
{
float X = p.x;
float Y = p.y;
// if(abs(Y) < 0.0001f) return PI_2 * sign(X);
float R = sqrt(X * X + Y * Y);
float sigma = asin(X / R);
// REFLECT IF η<0: ensure π/2<|σ|<π
sigma = (Y < 0.0f) ? ((sigma >= 0.0f) ? (PI - sigma) : -(PI + sigma)) : sigma;
return sigma;
}
// Maxwell-flow time at a point p=(X,Y) outside the unit circle.
// Evaluated using elliptic functions and Carlson form iterations.
// Ill-conditioned near X-axis.
float mxwTimeAt(float2 p)
{
float X = p.x;
float Y = abs(p.y); // MirrorY: We're in the upper court, kids.
float X2 = X * X;
float Y2 = Y * Y;
float R2 = X2 + Y2;
float R = sqrt(R2);
float sigma = asin(X / R);
sigma = mxwClampSigmaPosY(sigma); // Safety first.
float etaInf = Y - Y / R2;
float etaInf2 = etaInf * etaInf;
float etaSq4 = etaInf2 + 4.f;
float m = 4.f / etaSq4;
float2 FE = EllipticFE(sigma, m);
float F = FE.x;
float E = FE.y;
float x = (etaSq4 * E - (etaSq4 - 2.f) * F) / (2.f * sqrt(etaSq4));
float t = x - X;
return t;
}
// Maxwell-flow time at a point given by (σ,η∞).
// Ill-conditioned near X-axis.
// Assumption: |σ|<π/2 (Y>0)
float mxwTimeAtSigmaEta(float sigma, float eta)
{
sigma = mxwClampSigmaPosY(sigma); // Safety first.
eta = abs(eta); // MirrorY: We're in the upper court, kids.
float eta2 = eta * eta;
float etaSq4 = eta2 + 4.f;
float m = 4.f / etaSq4;
float2 FE = EllipticFE(sigma, m);
float F = FE.x;
float E = FE.y;
float x = (etaSq4 * E - (etaSq4 - 2.f) * F) / (2.f * sqrt(etaSq4));
float cosS = cos(sigma);
float Y = 0.5f * (eta + sqrt(eta2 + 4.f * cosS * cosS));
float X = Y * tan(sigma); // Ill-cond. near X axis: (Y→0)*(tan(σ→−π/2)→−∞)
float t = x - X;
return t;
}
// Maxwell-flow position given (σ,η∞).
// Badly behaved near X axis.
// Inputs:
// sigma - Angle σ((X,Y))∈(π,-π), where σ=0 corresponds to the +y axis.
// eta - Solution curve value, η_∞. Signed η values allowed.
// Returns: Position p=(X,Y) outside the unit circle.
float2 mxwPositionAtSigmaEta(float sigma, float eta)
{
float etaAbs = abs(eta);
float eta2 = eta * eta;
float cosS = cos(sigma);
float Y = 0.5f * (etaAbs + sqrt(eta2 + 4.f * cosS * cosS)); // Y>0 for abs(eta)
float sgnEta = sign(eta);
float X = sgnEta * Y * tan(sigma); // Ill-cond. near X axis: (Y→0)*(tan(σ→−π/2)→−∞)
Y *= sgnEta;
return float2(X, Y);
}
// Maxwell arctan monotone preconditioner, M, and its derivative, D, at (t, η∞).
// Returns float2(M,D).
float2 mxwMD_arctan(float t, float eta)
{
// REUSE S(eta):
float eta2 = eta * eta;
float S = 0.5f * (eta + (eta2 + 2.f) / sqrt(eta2 + 4.f)); // Use abs(eta) to make it symmetric top/bottom
float twoOverPI = 0.636619772367581343f;
float M = +twoOverPI * atan(-t / S);
float D = -twoOverPI / (S + t * t / S);
return float2(M, D);
}
// Maxwell tanha monotone preconditioner, M, and its derivative, D, at (t, η∞).
// Returns float2(M,D).
float2 mxwMD_tanha(float t, float eta)
{
// REUSE S(eta):
float eta2 = eta * eta;
float S = 0.5f * (eta + (eta2 + 2.f) / sqrt(eta2 + 4.f)); // Use abs(eta) to make it symmetric top/bottom
// Optimized magic constants
float a = 0.956764f;
float s = 1.05782f;
float c = 1.f / (PI_2 * s * S); // scaling (slope adjustment)
float T = -c * t; // Scaled time
float absT = abs(T);
float absTa = pow(absT, a);
float tanha = sign(T) * tanh(absTa);
float M = tanha;
// DERIV:
float coshTa = cosh(absTa);
float D = -a * c * absTa / (absT + 0.0001) / (coshTa * coshTa); // Regularize 1/|t| singularity
return float2(M, D);
}
// Maxwell monotone preconditioner's magic α* value (for η≥0).
float mxwMD_alphaStar(float eta)
{
float n = eta; // ouch, lazy
float n2 = n * n;
float squirt = sqrt(n2 + 4.f);
float alpha = 2.f * n * squirt / (n2 + n * squirt + 2.f);
return alpha;
}
// Maxwell monotone preconditioner M, and its derivative, D, at specified (t, η∞).
// Returns float2(M,D).
float2 mxwMD(float t, float eta)
{
float absEta = abs(eta);
float2 arctanMD = mxwMD_arctan(t, absEta);
float2 tanhaMD = mxwMD_tanha(t, absEta);
float alpha = mxwMD_alphaStar(absEta);
float2 MD = alpha * arctanMD + (1.f - alpha) * tanhaMD;
return MD;
}
// The notorious t'(σ)=dt/dσ integrand of Maxwell's integral for t(σ).
float mxwTimeDerivSigmaEta(float sigma, float eta)
{
float n = eta; // reuse lazily
float n2 = n * n;
float cosS = cos(sigma);
float cos2 = cosS * cosS;
float squirt = sqrt(n2 + 4.f * cos2);
return -0.5f * (n + squirt) * (n + squirt) / (cos2 * squirt);
}
/////////////////////////////////////////////////////////////////////////
/// Given Maxwell time tBar and η, solve for σ s.t. t(σ,η)=tBar. ///
/// Uses Newton's method with monotonic preconditioning, to solve ///
/// M(t(σ,η))=M(tBar) ⇔ f(σ)=fBar ///
/////////////////////////////////////////////////////////////////////////
// NOTE: Works for all tBar and eta signs.
// Returns |σ|<π/2 for η>0.
// Returns π/2<|σ|<π for η<0.
float mxwSigmaSolver(float tBar, float eta)
{
float etaAbs = abs(eta);
// EVAL fBar:
float2 Fbar = mxwMD(tBar, etaAbs); // sloth: don't need derivative, D
float fBar = Fbar.x; // fBar = M(tBar); Range: (-1,1)
////////////////////////
// NEWTON'S METHOD: //
////////////////////////
// INITIAL GUESS (straight line approximation)
float sigma = fBar * PI_2; // <π/2
sigma = mxwClampSigmaPosY(sigma); // Safety first
for (int iter = 0; iter < 20; iter++)
{
float t = mxwTimeAtSigmaEta(sigma, etaAbs);
float2 MD = mxwMD(t, etaAbs); // NOTE: eta>0 enforced
float f = MD.x;
float fDeriv = MD.y * mxwTimeDerivSigmaEta(sigma, etaAbs);
float dSigma = -(f - fBar) / fDeriv;
sigma += dSigma;
// ADAPTIVE NEWTON TOL:
float SIGMA_NEWTON_TOL = 0.0000002f; // Newton solver tol for σ(t)
// float tolExp = 5.0f*smoothstep(0.01f, 0.5f, etaAbs);
// SIGMA_NEWTON_TOL *= pow(2.f, tolExp);
// if(etaAbs > 1.f) SIGMA_NEWTON_TOL = 0.00001f;
// SIGMA_NEWTON_TOL = 0000002f; // DEBUG: OVERRIDE ADAPTIVE TOL
// const float SIGMA_NEWTON_TOL = 0.000001f; // Newton solver tol for σ(t)
if (abs(dSigma) < PI_2 * SIGMA_NEWTON_TOL)
break;
}
// REFLECT IF η<0: ensure π/2<|σ|<π
if (eta < 0.f)
{
if (sigma >= 0.f)
sigma = PI - sigma;
else
sigma = -(PI + sigma);
}
return sigma;
}
// Raw monotone-preconditioned Newton solve for position at (t, η).
// Poor near η≅0.
float2 mxwPositionSolveAtTimeEta(float t, float eta)
{
eta = (eta != 0.0f) ? eta : 0.0001f; // Please stay off the line folks.
if (abs(eta) < 0.0001f)
eta = 0.0001f * sign(eta);
float sigma = mxwSigmaSolver(t, eta); // RETURNS σ∈(π,-π). Bad for large |t|.
return mxwPositionAtSigmaEta(sigma, eta); // Poor near η≅0.
}
// Raw monotone-preconditioned Newton solver for position advected for dt time.
float2 mxwAdvectPointRaw(float2 p, float dt)
{
float eta = mxwEtaInf(p);
float t = mxwTimeAt(p); // BUG: Huge values near eta~0
float tNew = t + dt;
return mxwPositionSolveAtTimeEta(tNew, eta); // BUG: Huge values of tNew near η~0.
}
// Advects point through unit-cylinder flow but transitions to "pure translation"
// outside a maximum T/Dist range to avoid distortion near the x-axis.
// Parameters:
// p - Starting position.
// dt - Advection duration.
// tCutoff - A positive time/distance, e.g., 20, beyond which to use the translation approximation.
// Returns:
// Final advected position.
float2 mxwAdvectPointTCutoff(float2 p, float dt, float tCutoff)
{
float eta = mxwEtaInf(p);
float t = mxwTimeAt(p); // BUG: Huge/Poor values near eta~0
float tNew = t + dt;
float2 pNew;
if (abs(tNew) < tCutoff)
{
pNew = mxwPositionSolveAtTimeEta(tNew, eta); // BUG: Huge values of tNew near η∞~0.
}
else
{
float tSafe = sign(tNew) * tCutoff;
float tSlide = tNew - tSafe;
pNew = mxwPositionSolveAtTimeEta(tSafe, eta);
pNew.x -= tSlide;
}
return pNew;
}
// Advects point for time dt through flow around unit circle.
float2 mxwAdvectPointUnpatched(float2 p, float dt)
{
return mxwAdvectPointTCutoff(p, dt, 20.f); // Hard-coded TCutoff value
}
// Advects p for dt time using a patched version to handle singularities near η~0.
float2 mxwAdvectPointUnitPatched(float2 p, float dt, bool applyPatch)
{
float etaInf_p = mxwEtaInf(p); // Evaluate η∞ after tine insertion.
float2 pNew = p;
float etaMax = 0.05f; // Patch size threshold
if (applyPatch && abs(etaInf_p) < etaMax)
{
// EXTRAPOLATE INWARD ON y=0 PARAMETERIZATION'S SEAM:
bool inWedge = (acos(abs(normalize(p).x)) < 0.05f);
float xCrease = sign(dt) * (0.5f + abs(dt));
bool inPatchRegion =
(inWedge && ((dt < 0.f && (p.x >= +1.f || p.x < xCrease)) || (dt > 0.f && (p.x <= -1.f || p.x > xCrease))));
if (inPatchRegion)
{
// Extrapolate from off-line y∈[-H,H]
float H = etaMax;
float dy = (p.y < 0.f ? -H : H);
float2 pPerturb = float2(p.x, dy);
float2 pPerturbNew = mxwAdvectPointUnpatched(pPerturb, dt); // [1st advection eval]
float2 u = pPerturbNew - pPerturb;
// u.y SINGULARITY CORRECTION: Simple and fast linear model
u.y *= abs(p.y) / H;
// Optional: cubic correction for u.y (meh):
if (false)
{
float uy0 = u.y;
float uy1 = 0.f; // Placeholder for secondary value if needed
float a = (8.f * uy0 - uy1) / (6.f * H);
float b = (uy1 - 2.f * uy0) / (6.f * H * H * H);
float y = abs(p.y);
u.y = a * y + b * y * y * y;
}
// Fancy u.x curvature correction:
if (false)
{
float2 pPertur2 = float2(p.x, dy * 2.f);
float2 pPertur2New = mxwAdvectPointUnpatched(pPertur2, dt); // [2nd advection eval]
float2 u2 = pPertur2New - pPertur2;
// Parabolic model for u.x(y) correction
float ux1 = u.x;
float ux2 = u2.x;
float a = (ux2 - ux1) / (3.f * H * H);
float b = (4.f * ux1 - ux2) / 3.f;
float y = p.y;
u.x += (a * y * y + b - u.x);
}
pNew = p + u;
return pNew;
}
}
pNew = mxwAdvectPointUnpatched(p, dt); // Default unpatched flow
return pNew;
}
// Maps point p to account for 2D incompressible tine insertion.
float2 mxwApplyTineInsertionMap(float2 p, float r)
{
if (dot(p, p) == 0.f)
p = 0.000001f * float2(1.f, 1.f); // avoids R=0 and 0/0 result.
float R = length(p); // current radius
float S = sqrt(R * R + r * r);
return (S / R) * p;
}
// Maps point p to account for 2D incompressible tine removal.
float2 mxwApplyTineRemovalMap(float2 p, float r)
{
float R = length(p); // current radius
float S = sqrt(max(0.f, R * R - r * r)); // avoid spurious NaNs
return (S / R) * p;
}
// Advects point p around a unit-cylinder flow considering tine insertion/removal and singularity patch.
float2 mxwAdvectPointUnit(float2 p, float dt, bool insertRemoveTine, bool applyPatch)
{
float2 p0 = p; // initial position
// SETUP CUSP PATCH:
bool cuspPatch = false;
float cuspH = min(0.02f, 0.005f * abs(dt));
if (insertRemoveTine && applyPatch)
{
if (abs(p0.y) < cuspH && (sign(dt) * p0.x < 0.f || sign(dt) * p0.x > abs(dt) * 0.8f))
{
float dy = (p0.y < 0.f ? -cuspH : cuspH);
p.y = dy; // perturb compute location
cuspPatch = true; // apply blend at end to u.y
}
}
float2 pPreCusp = p; // for computing displacement, u
// Handle tine insertion: inflate coordinates for R=1+epsRegularization circle.
if (insertRemoveTine)
{
float epsR = (applyPatch ? 0.0002f : 0.f); // pull out slightly more paint during tine removal
p = mxwApplyTineInsertionMap(p, 1.f + epsR);
}
// Advection around patched unit cylinder.
// float2 pNew = mxwAdvectPointUnitPatched(p, dt, applyPatch);
float2 pNew = mxwAdvectPointUnpatched(p, dt); // DEBUG/NEW (only use cuspPatch if applyPatch)
// Handle tine removal: deflate coordinates.
if (insertRemoveTine)
{
pNew = mxwApplyTineRemovalMap(pNew, 1.f);
}
// PATCH CUSP:
if (cuspPatch)
{
float2 u = pNew - pPreCusp;
u.y *= abs(p0.y) / cuspH;
pNew = p0 + u; // filtered u.y
}
return pNew;
}
// Advects point p around a cylinder of radius R for flow distance L in -x direction.
float2 mxwAdvectPointScaled(float2 p, float L, float R, bool insertRemoveTine, bool applyPatch)
{
// Scale to unit-radius world.
p /= R;
float dt = L / R; // distance (time) in unit-radius world.
// Advect point.
float2 pNew = mxwAdvectPointUnit(p, dt, insertRemoveTine, applyPatch);
// Return unscaled point.
return R * pNew;
}
/******************************************************************************/
/*** [END] MAXWELL FLOW PROBLEM AND MONOTONE-PRECONDITIONED NEWTON SOLVER ***/
/******************************************************************************/
/*****************************************************************************/
/*** [BEGIN] REVERSE-DRIFT FIELD (RDF) IMPLEMENTATIONS **/
/*****************************************************************************/
// Line RDF: Matched series approximation to 1D drift.
// Input: height above flow line, and cylinder radius, r.
float rdLine1DS(float height, float r)
{ /// DEBUG/NEW
float y = height / r;
float dx = 0.f; // drift = 5.f * r / (0.2f + pow(eta, 3.f));
float eps = 0.002f;
y = sqrt(y * y + eps * eps); // y_eps
if (y > 6.50416858646776f)
{
float h = 1.f / y;
float h2 = h * h;
float h3 = h * h2;
float poly = ((((6615.f * h2 - 1960.f) * h2 + 600.f) * h2 - 192.f) * h2 + 64.f);
dx = -(PI / 256.f) * h3 * poly;
}
else if (y > 0.8220844420096408f)
{
float dy = y - 2.3f;
float num = (((((-1.01973e-7f * dy + 1.77539e-6f) * dy - 0.0068946f) * dy - 0.043355f) * dy - 0.0937878f) * dy - 0.0413839f);
float den = ((((((0.00875686f * dy + 0.115772f) * dy + 0.665761f) * dy + 2.08615f) * dy + 3.66495f) * dy + 3.26035f) * dy + 1.0f);
dx = num / den;
}
else
{
float z = y * y;
float A = ((((0.00015811568f * z - 0.00096477622f) * z + 0.0076619254f) * z - 0.16369764f) * z - 0.039720771f);
float B = ((((-0.00018775463f * z + 0.0010681152f) * z - 0.0073242188f) * z + 0.09375f) * z + 0.5f);
dx = A + B * log(y);
}
return 2.0f * r * dx; // UNDO UNIT-RADIUS XFORM & *2 for full drift.
}
// Returns rdLine1DS series implementation of 1D drift. Provided for "backwards convenience."
float drift1D(float height, float r)
{
return rdLine1DS(height, r);
}
// Reverse-drift field (RDF) for a line segment (a→b) computed with Mixwell solver.
// Tine insertion and removal is assumed, so the RDF can be evaluated everywhere.
// The Mixwell solver singularity is patched.
//
// Parameters:
// p : Position to evaluate RDF
// r : Tine radius
// a : Line-segment start
// b : Line-segment end
float2 rdSegmentM(float2 p, float r, float2 a, float2 b)
{
if (length(a - b) == 0.f)
return float2(0.f, 0.f);
// Reverse of a→b motion is a←b:
// So shift b to the origin (start):
p -= b;
float2 M = a - b; // Reverse motion vector
// Rotate to x-axis flow:
float L = length(M); // assume dt=L/r>0 always
float2 newXDir = normalize(M);
float cs = newXDir.x;
float sn = newXDir.y;
// float2x2 Q = (float2x2)(cs, -sn, sn, cs); // rotation of X to B
// float2x2 QT = (float2x2)(cs, sn, -sn, cs); // rotation of B to X
if (dot(p, p) == 0.)
return M; // point at cusp tip.
float2 p0 = p;
// Rotate to x-axis flow orientation //p = mul(Q, p)
p = float2(cs * p.x + sn * p.y, -sn * p.x + cs * p.y);
bool insertRemoveTine = true;
bool applyPatch = true;
float2 pNew = mxwAdvectPointScaled(p, L, r, insertRemoveTine, applyPatch);
// Unrotate back to B flow orientation //pNew = mul(QT, pNew)
pNew = float2(cs * pNew.x - sn * pNew.y, +sn * pNew.x + cs * pNew.y);
return (pNew - p0) + M;
}
// 1D reverse drift for a line (e.g., x-axis) computed using Elliptic integrals (E).
// The singularity is regularized.
//
// Parameters:
// height - Height above flow centerline
// r - Cylinder radius
float rdLine1DE(float height, float r)
{
// Convert to unit-radius problem:
float y = height / r; // = η∞
float y2 = y * y; // = η∞²
// Rectify y and regularize singularity
float eps = 0.002f;
y2 += eps * eps; // y2=y²+ε²
y = sqrt(y2); // y≥ε
// Complete elliptic integrals: K[m]=F(π/2,m) and E[m]=E(π/2,m)
float m = -4.0f / y2;
float2 FE = EllipticFE(PI_2, m);
float K = FE.x; // K[m]=F(π/2,m)
float E = FE.y; // E[m]=E(π/2,m)
// 1D drift
float dx = (y2 * E - (y2 + 2.f) * K) / y;
return r * dx; // Scale to r-radius problem
}
// vec2 rdLineE() {}
// cylFlowVelWorld: Velocity at point p in World frame for unit cylinder flow.
// Assumes unit radius cylinder at origin; unit translation speed along x; and v=(0,0) at ∞.
// Input: float2 p - 2D point.
// Output: float2 v - Velocity vector in World frame.
float2 cylFlowVelWorld(float2 p)
{
p = projectOutsideUnitDisk(p); // Safety first! (for timestepper)
float X2 = p.x * p.x; // Precompute X²
float Y2 = p.y * p.y; // Precompute Y²
float R2 = X2 + Y2; // R² = X² + Y²
float invR4 = 1.0f / (R2 * R2); // Inverse R⁴
return float2((X2 - Y2) * invR4, 2.f * p.x * p.y * invR4);
}
// cylFlowVelBody: Velocity at point p in Body frame for unit cylinder flow.
// Assumes unit radius cylinder at origin; unit translation speed along x; and v=(-1,0) at ∞.
// Input: float2 p - 2D point. Output: float2 - Velocity vector in Body frame.
float2 cylFlowVelBody(float2 p)
{
return cylFlowVelWorld(p) - float2(1.f, 0.f);
}
// cylFlowVelWorldArbDir: Velocity at point p in World frame for unit cylinder flow.
// Assumes unit radius cylinder at origin; unit translation speed along arbitrary direction vc; and v=(0,0) at ∞.
// Input: float2 p - 2D point.
// float2 vc - Unit velocity of cylinder.
// Output: float2 - Velocity vector in World frame.
float2 cylFlowVelWorldArbDir(float2 p, float2 vc)
{
float2 d = normalize(vc); // unit direction of cylinder motion
float2 pRef = float2(d.x * p.x + d.y * p.y, -d.y * p.x + d.x * p.y); // Rotate p to reference frame
float2 vRef = cylFlowVelWorld(pRef); // Velocity in reference frame
return float2(d.x * vRef.x - d.y * vRef.y, d.y * vRef.x + d.x * vRef.y); // Rotate velocity back to World frame
}
// Verification: Brute-force time-stepped (TS) advection approximation to 1D drift. (Body frame calculation)
float rdLine1DTBody(float height, float r)
{
float y = height / r; // UNIT-RADIUS TRANSFORM
// RELEASE PARTICLE UPSTREAM (x=Xmax) and ADVECT UNTIL x<-Xmax
// COMPARE DRIFT TO PURE TRANSLATING PARTICLE:
float Xmax = 1000.f; // big enough to avoid negative drift
float2 p0 = float2(Xmax, y); // Particle init at (0,y) --> (pureTranslation-dx, ~y)
float2 p = p0;
int i = 0;
float time = 0.0f;
while (p.x > -Xmax && i < 3000)
{
float dt = 0.01f * length(p); // adaptive stepsize
float2 vp = cylFlowVelBody(p); // initial velocity; (-1,0) at ∞
float2 pmid = p + vp * dt * 0.5f; // midpoint position
float2 vmid = cylFlowVelBody(pmid); // midpoint velocity
p += vmid * dt; // midpoint method step
time += dt;
i++;
}
float xPureTranslation = Xmax - time;
float dx = -(p.x - xPureTranslation); // 1D drift
return dx * r; // UNDO UNIT-RADIUS TRANSFORM
}
// Verification: Brute-force time-stepped advection approximation to 1D drift. (World frame calculation)
float rdLine1DTWorld(float height, float r)
{
float y = height / r; // UNIT-RADIUS TRANSFORM
// ADVECT PARTICLE at (0,y) AROUND CYLINDER MOVING FROM (x=-Xmax to +Xmax).
// X DISPLACEMENT IS DRIFT.
float Xmax = 1000.f; // big enough to avoid negative drift
float2 p0 = float2(0.f, y); // Particle init at (0,y) --> (-dx, ~y)
float2 c = float2(-Xmax, 0.f); // Cylinder position
float2 vc = float2(1.f, 0.f); // Cylinder velocity (constant for 1D drift calc)
float2 p = p0;
int i = 0;
float time = 0.0f;
while (c.x < Xmax && i < 3000)
{
float2 q = p - c; // relative particle position
float dt = 0.01f * length(q); // adaptive stepsize
float2 vp = cylFlowVelWorld(q); // initial velocity; (-1,0) at ∞
float2 pmid = p + vp * dt * 0.5f; // midpoint position
float2 cmid = c + vc * dt * 0.5f; // midpoint cylinder position
float2 qmid = pmid - cmid; // midpoint relative particle position
float2 vmid = cylFlowVelWorld(qmid); // midpoint velocity
p += vmid * dt; // midpoint method step
c += vc * dt; // advance cylinder
time += dt;
i++;
}
float dx = -p.x; // 1D drift
return dx * r; // UNDO UNIT-RADIUS TRANSFORM
}
// RDF at p for a circle(r) passing along an infinite line in (dir)ection passing through origin.
float2 rdLine(float2 p, float r, float2 dir)
{
dir = normalize(dir);
float2 n = rot90(dir); // float2(-dir.y, +dir.x);
float height = dot(p, n);
float2 rdf = rdLine1DS(height, r) * dir;
return rdf;
}
// Min distance from y to a tine of a comb with pitch hy, rooted at the origin.
float minCombDist(float y, float hy)
{
return hy * abs(frac(y / hy - 0.5f) - 0.5f);
}
// RDF for Nonpareil pattern.
// Inputs:
// p : Evaluation point.
// r : Tine radius.
// dir : Combing direction.
// combGap : Spacing between tines.
// nPasses : Number of nearby tines to approximate drift by.
float2 rdNonpareil(float2 p, float combR, float2 combDir, float combGap, int nPasses)
{
float2 dir = normalize(combDir);
float2 n = rot90(dir);
float y = dot(p, n);
float passGap = float(nPasses) * combGap; // wider spacing between passes for nicer falloff + overlap
float2 rdf = float2(0.f, 0.f);
for (int k = 0; k < nPasses; k++)
{
float yOffset = float(k) * combGap;
float height = minCombDist(y + yOffset, passGap);
rdf += rdLine1DS(height, combR) * dir;
}
return rdf;
}
// RDF for Nonpareil pattern with comb gap noise.
// Inputs:
// p : Evaluation point.
// r : Tine radius.
// dir : Combing direction.
// combGap : Spacing between tines.
// combGapNoise : Amplitude of tine position variation (± this value).
float2 rdNonpareilNoisy(float2 p, float r, float2 dir, float combGap, float combGapNoise)
{
dir = normalize(dir);
float2 n = rot90(dir);
float y = dot(p, n);
float dyMax = r * RDSEGMENT_MASK_R_FACTOR + combGap;
int dkMax = (int) (ceil(dyMax / combGap));
int k0 = (int) (round(y / combGap)); // closest (unpeturbed) tine.
float2 rdf = float2(0.f, 0.f);
for (int k = k0 - dkMax; k <= k0 + dkMax; k++)
{
float noisek = 2.f * iqhash(float(k)) - 1.f; // ∈[-1,1)
float yk = combGap * float(k) + combGapNoise * noisek;
////if(iqhash(noisek*FUDGE) > 0.12) // Optional tine dropout (skip for now)
rdf += rdLine1DS(abs(y - yk), r) * dir;
}
return rdf;
}
// RDF for Gel-Git pattern.
float2 rdGelGit(float2 p, float combR, float2 combDir, float combGap, int nPasses)
{
float2 dir = normalize(combDir);
float2 n = rot90(dir);
float y = dot(p, n);
float passGap = float(nPasses) * combGap * 2.f; // wider spacing between passes for nicer falloff + overlap
float2 rdf = float2(0.f, 0.f);
// GIT ("left")
for (int k = 0; k < nPasses; k++)
{
float yOffset = (0.5f + float(k)) * combGap * 2.f;
float height = minCombDist(y + yOffset, passGap);
rdf -= rdLine1DS(height, combR) * dir;
}
// GEL ("right")
for (int k2 = 0; k2 < nPasses; k2++)
{
float yOffset = float(k2) * combGap * 2.f;
float height = minCombDist(y + yOffset, passGap);
rdf += rdLine1DS(height, combR) * dir;
}
return rdf;
}
// Returns A*sin(K*X+P) given W=(A,K,P).
float waveEval(float X, float3 W)
{
return W.x * sin(W.y * X + W.z);
}
// Returns unit vector tangent at X for f(X)=A*sin(K*X+P) given W=(A,K,P).
float2 waveUnitTangent(float X, float3 W)
{
float Df = W.x * W.y * cos(W.y * X + W.z); // f'(X)
float2 T = float2(1.f, Df);
return normalize(T);
}
// Time-stepped advection approximation to 2D drift from a sine wave motion. (World frame calculation)
// Sine wave: y0 + amp*sin(K x + phase)
// Input: float3 W=(A,K,P) describes the partial wave A*sin(K*X+P)
// rightComb: true passes right (+x), false passes left (-x) direction.
float2 rdSineT(float2 p, float r, float waveY0, float3 W, bool rightComb)
{
float waveAmp = W.x;
float waveK = W.y;
float wavePhase = W.z;
p.y -= waveY0; // TRANSLATE TO Y=0 CENTERLINE
// UNIT-RADIUS TRANSFORM:
p /= r;
waveAmp /= r;
waveK *= r; // K'*x'=K*x invariance with x'=x/r implies that K'=K*r
// TRANSLATE PARTICLE TO X=0 (move offset into phase):
wavePhase += waveK * p.x;
p.x = 0.f; // now X offset
float y0 = p.y; // INIT POS: p=(0,y0)
// WAVE PARAMS:
W = float3(waveAmp, waveK, wavePhase); // scaled wave params
float rdir = (rightComb ? -1.f : +1.f); // reverse x-dir for combing.
// ADVECT PARTICLE at p AS CYLINDER MOVES ON SINEWAVE FROM (x=-Xmax to +Xmax).
float Xmax = 50.f; // big enough to avoid negative drift
float2 p0 = float2(0.f, y0); // Particle init at (0,y0) --> p
float2 c = float2(-Xmax * rdir, 0.f); // Cylinder position (bogus y value)
float2 vc; // Cylinder velocity (unit vector)
int i = 0;
float time = 0.f;
while ((rdir * c.x < Xmax) && i < 5000)
{
c.y = waveEval(c.x, W); // correct cylinder y (project onto sine)
vc = waveUnitTangent(c.x, W) * rdir; // Unit velocity vector (vc.x > 0)
float2 q = p - c; // relative particle position
float dt = 0.05f; // HQ:0.01 SQ:0.05
float2 vp = cylFlowVelWorldArbDir(q, vc); // particle velocity
float2 pmid = p + vp * dt * 0.5f; // midpoint position
float2 cmid = c + vc * dt * 0.5f; // midpoint cylinder position
float2 vcmid = waveUnitTangent(cmid.x, W) * rdir;
float2 qmid = pmid - cmid; // midpoint relative particle position
float2 vmid = cylFlowVelWorldArbDir(qmid, vcmid); // midpoint velocity
// MIDPOINT METHOD:
p += vmid * dt;
c += vcmid * dt;
p = projectOutsideUnitDisk(p - c) + c; // Avoid losing particles inside disk
time += dt;
i++;
}
float2 rdrift = p - p0; // reverse drift
return rdrift * r; // UNDO UNIT-RADIUS TRANSFORM
}
// DEBUG TEST: 1D drift computed using degenerate sine wave 2D drift. <VERIFIED!>
float debugRDrift1D_sine(float height, float r)
{
float2 rdrift = rdSineT(float2(0.f, height), r, 0.f, float3(0.f, 1.f, 0.f), true);
return rdrift.x;
}
// Combs a→b (performs insert/remove flows at begin/end)
// Time-stepped implementation (T).
float2 rdSegmentT(float2 p, float r, float2 a, float2 b)
{
// UNIT-RADIUS TRANSFORM:
float invR = 1.f / r;
p *= invR;
a *= invR;
b *= invR;
// REVERSE COMB SETUP: b→a
float2 rdir = a - b;
float L = length(rdir); // units of r
float dist2Segment = sdSegment(p, a, b);
// ADVECT PARTICLE at p AS CYLINDER MOVES FROM b to a.
float2 p0 = p; // Particle init (USED TO GET RevDrift: p-p0)
float2 c = b; // Cylinder init
float2 vc = normalize(rdir); // Cylinder unit velocity
int i = 0;
float time = 0.f;
// IF REMOVE TINE, INFLATE COORDS:
p = b + mxwApplyTineInsertionMap(p - b, 1.f); // WORLD FRAME
// Reverse-advect p as cylinder moves from b→a:
float tLeft = L;
while (tLeft > 0.000f)
{
float2 q = p - c; // relative particle position
float dt = 0.2f; // adaptive stepsize (near segment)
dt *= max(0.05f, dist2Segment);
if (dt > tLeft)
dt = tLeft;
float2 vp = cylFlowVelWorldArbDir(q, vc); // particle velocity
float2 pmid = p + vp * dt * 0.5f; // midpoint position
float2 cmid = c + vc * dt * 0.5f; // midpoint cylinder position
float2 qmid = pmid - cmid; // midpoint relative particle position
float2 vmid = cylFlowVelWorldArbDir(qmid, vc); // midpt velocity
// MIDPOINT METHOD:
p += vmid * dt; // position update
c += vc * dt; // advance cylinder
p = projectOutsideUnitDisk(p - c) + c; // Avoid losing particles inside disk
time += dt;
tLeft -= dt;
i++;
}
p = a + mxwApplyTineRemovalMap(p - a, 1.f); // WORLD FRAME; unit tine
float2 rdrift = p - p0; // reverse drift
return rdrift * r; // UNDO UNIT-RADIUS TRANSFORM
}
// Mixwell brush matrix for relative position, v, and radius eps. // NEW
float2x2 getMixwellBrushMatrix(float2 v, float eps)
{
float R2 = dot(v, v);
float e2 = eps * eps;
float s = R2 + e2;
// invD = 1/(s^(3/2)) = 1/(s*sqrt(s))
float invSqrtS = rsqrt(s);
float invS = invSqrtS * invSqrtS;
float invD = invS * invSqrtS;
float invR = rsqrt(max(R2, 1e-20f)); // guard 1/R at v≈0
float R = R2 * invR; // sqrt(R2); WIN: Saves an XU instruction!
// Afd = 1 - R*(R2 + 2e2)/d == 1 - sqrt(R2)*(R2 + 2e2)*invD
float Afd = 1.0f - R * (R2 + 2.0f * e2) * invD;
// Bfd = e2/(R*d) == e2*(1/R)*invD
float Bfd = e2 * invR * invD;
// K = Afd*I + Bfd*(v v^T)
float xx = Bfd * v.x * v.x;
float xy = Bfd * v.x * v.y;
float yy = Bfd * v.y * v.y;
return float2x2(Afd + xx, xy, xy, Afd + yy);
}
// Reverse drift through Mixwell flow of radius-eps brush moving a→b. Computed using adaptive midpoint integration. //
// NEW
float2 rdMixwellBrushAdaptiveMidpoint(float2 p0, float eps, float2 a, float2 b)
{
// Step brush position backwards from b to a, advecting p to compute rev-drift (p-p0):
float2 uBrush = a - b; // rev brush step
float L = length(uBrush);
if (L <= 1e-20f)
return float2(0, 0);
float2 dir = uBrush / L; // rev brush dir
float distLeft = L; // brush motion left
float2 p = p0; // reverse-drifted point position (init: p0)
while (distLeft > 0.f)
{
float2 v = p - b; // rel pos from end brush goal
float R = length(v);
float dL = min(distLeft, 0.10f * max(eps, R)); // adaptive stepsize (spatially continuous errors)
distLeft -= dL;
// MIDPOINT STEP: dp = K(vmid) db
float2 db = dL * dir; // Δbrush
float2x2 K = getMixwellBrushMatrix(v, eps); // K(v)
float2 dp = mul(K, db); // Euler approx, dp = K(v) db
float2 vmid = v + 0.5f * (dp - db); // v' = p' - b' = v + 0.5*(dp - db)
float2x2 Kmid = getMixwellBrushMatrix(vmid, eps); // K(vmid)
float2 dpmid = mul(Kmid, db); // Midpoint step
p += dpmid; // update point (rd midpoint approx)
b += db; // update brush location
}
return p - p0; // reverse-drift displacement
}
// Reverse drift for line segment (a→b) using a time-stepped Mixwell Brush (B) of radius r. // NEW
float2 rdSegmentMBrush(in float2 p, in float r, in float2 a, in float2 b)
{
return rdMixwellBrushAdaptiveMidpoint(p, r, a, b);
}
// Combs a→b (performs insert/remove flows at begin/end)
float2 rdSegment(float2 p, float r, float2 a, float2 b)
{
// Optional: LOCAL MASK (zero beyond dMax ∝ r):
float d = sdSegment(p, a, b);
float dMax = RDSEGMENT_MASK_R_FACTOR * r; // ADJUST RDSEGMENT_MASK_R_FACTOR TO TASTE
if (d >= dMax)
return float2(0.f, 0.f);
float mask = 1.f - smoothstep(0.5f * dMax, dMax, d); // blend-to-zero factor.
// return mask * rdSegmentM(p, r, a, b); // Mixwell Newton solver
return mask * rdSegmentMBrush(p, r, a, b); // Mixwell Brush solver // NEW
}
// Combs a→b (performs insert/remove flows at begin/end), with reverse option.
// @param reverse - Swaps a↔b.
float2 rdSegmentRev(float2 p, float r, float2 a, float2 b, bool reverse)
{
float2 head = (reverse ? b : a);
float2 tail = (reverse ? a : b);
return rdSegment(p, r, head, tail);
}
float2 rdTriangle(float2 p, float r, float2 p0, float2 p1, float2 p2)
{
float2 pInit = p; // copy to compute rdrift
float2 P[4] = { p0, p1, p2, p0 };
float sdf = sdTriangle(p, p0, p1, p2);
float sdfMax = 5.f * r;
if (sdf > sdfMax)
return p - pInit;
// Masking factor to avoid discontinuity:
float masking = 1.f - smoothstep(sdfMax / 2.f, sdfMax, abs(sdf));
// JAGGED CASE (as rdrift superposition - symmetric but nonphysical)
float2 u = float2(0.f, 0.f);
u += rdSegmentRev(p, r, P[0], P[1], false);
u += rdSegmentRev(p, r, P[1], P[2], false);
u += rdSegmentRev(p, r, P[2], P[0], false);
p += u;
return (p - pInit) * masking;
}
// FOR ANALYSIS: Segment kernel (r=1) for L-length brush stroke from a=(0,0) to b=(L,0).
float2 rdSegmentKernel(float2 p, float L)
{
float2 a = float2(0.f, 0.f);
float2 b = float2(L, 0.f);
float r = 1.f;
// REVERSE COMB SETUP: b→a
float2 rdir = a - b;
// ADVECT PARTICLE at p AS CYLINDER MOVES FROM b to a.
float2 p0 = p; // Particle init (USED TO GET RevDrift: p-p0)
float2 c = b; // Cylinder init
float2 vc = normalize(rdir); // Cylinder unit velocity
// IF REMOVE TINE, INFLATE COORDS:
p = b + mxwApplyTineInsertionMap(p - b, 1.f); // WORLD FRAME
// Reverse-advect p as cylinder moves from b→a:
float2 q = p - c; // relative particle position
float dt = L;
float2 vp = cylFlowVelWorldArbDir(q, vc); // particle velocity
p += vp * dt; // position update
c += vc * dt; // advance cylinder
p = a + mxwApplyTineRemovalMap(p - a, 1.f); // unit tine
float2 rdrift = p - p0; // reverse drift
return rdrift;
}
/**
* Triangle Wave RDF.
*
* Uses rdSegments to construct a triangle wave parameterized like a sine wave
* at the origin with user specified direction and parameters.
* Uses sparse evaluation of masked rdSegments only near p.
*
* Inputs:
* p - Evaluation point.
* r - Tine radius.
* dir - Direction of wave centerline.
* L - Wavelength
* A - Signed amplitude. Use A=0 for broken and dashed lines.
* broken - Broken segments using doubly reversed rdSegments.
* dashed - Draws only every other segment.
*/
float2 rdTriWave(float2 p, float r, float2 dir, float L, float A, bool broken, bool dashed)
{
float2 p0 = p; // Initial position for reverse drift calc (p - p0)
dir = normalize(broken ? -dir : dir); // Change direction for reversed rdSegments flow trend
// Project p onto line --> c and get coords for p as (x, y)
float2 c = dot(p, dir) * dir;
float x = dot(c, dir);
float y = length(p - c); // Unsigned distance to line
float segR = sqrt(A * A + L * L / 16.0f);
float R = r * RDSEGMENT_MASK_R_FACTOR + segR; // Bounding radius of point w.r.t. qi
float H = L * 0.5f;
// Guard against a point outside the influence band (R*R - y*y < 0 -> sqrt is NaN) and a
// degenerate wavelength (H <= 0). Casting NaN/Inf to int for the loop bounds below is
// undefined behavior; under WebGPU/Dawn it produces a runaway loop that hangs the GPU
// (DXGI_ERROR_DEVICE_HUNG). Outside the band no segment contributes, so return zero drift.
float disc = R * R - y * y;
if (disc <= 0.0f || H <= 0.0f)
return float2(0.0f, 0.0f);
float Rx = sqrt(disc); // Influence "radius" about x on line
int iR = (int) ceil((x + Rx) / H); // max i
int iL = (int) floor((x - Rx) / H); // min i
float2 vd = 0.5f * H * dir;
float2 vp = A * rot90(dir);
for (int i = iR; i >= iL; i--)
{ // Reverse sweep
float si = (abs(i) % 2 == 0) ? -1.0f : 1.0f; // -1 even | +1 odd
if (!dashed || (dashed && si < 0.0f))
{
float2 qi = dir * H * float(i);
float2 ai = qi - vd + si * vp;
float2 bi = qi + vd - si * vp;
p += rdSegmentRev(p, r, ai, bi, broken);
}
}
return p - p0; // Reverse drift displacement
}
float2 rdTriWaveComb(float2 p, float r, float2 dir, float L, float A, bool broken, bool dashed, float combGap)
{
dir = normalize(dir);
float2 n = rot90(dir);
float y = dot(p, n);
float R = A + r * RDSEGMENT_MASK_R_FACTOR; // 1D perp-influence radius of a wave centerline
if (combGap <= 0.0f)
return float2(0.0f, 0.0f); // degenerate spacing -> division below would be Inf (GPU hang)
int dk = (int) ceil(R / combGap); // Conservative range of wave indices influencing p
int kp = (int) round(y / combGap); // Closest wave index (where origin wave is k=0)
float2 rdf = float2(0.f, 0.f);
for (int k = kp - dk; k <= kp + dk; k++)
{ // Sum nearby(p) wave RDFs:
float yk = float(k) * combGap; // Perp offset of k'th centerline
float2 pk = p - yk * n; // Translate k'th wave to origin
rdf += rdTriWave(pk, r, dir, L, A, broken, dashed); // Accumulate RDF of k'th triwave
}
return rdf;
}
/// END/////////////////////////////////////////////////////////////////Part of Mixwell HLSL. Released under the MIT license.