Mixwell.cl
The full Mixwell library: brushes, particle drift, and Reverse-Drift Functions, along with the math helpers they depend on.
// SPDX-License-Identifier: MIT
// Copyright (c) 2026 Doug L. James and Ethan James
///////////////////////////////////////////////////////////////////////
//// Mixwell Library Functions /// @author Doug L James, 2026 ////////
///////////////////////////////////////////////////////////////////////
#define PI 3.14159265358979323846f // π
#define PI_2 1.57079632679489661923f // π/2
// Cutoff distance (as multiple of tine radius) beyond which rdSegment is zero.
__constant float RDSEGMENT_MASK_R_FACTOR = 10.0f;
__attribute__((overloadable)) float fract(float x) {
return x - floor(x);
}
__attribute__((overloadable)) float2 fract(float2 x) {
return (float2)(fract(x.x), fract(x.y));
}
__attribute__((overloadable)) float3 fract(float3 x) {
return (float3)(fract(x.x), fract(x.y), fract(x.z));
}
__attribute__((overloadable)) float4 fract(float4 x) {
return (float4)(fract(x.x), fract(x.y), fract(x.z), fract(x.w));
}
// Euclidean mod // NEW
float emod (float x, float n){ return fract(x / n) * n;}
float2 emod2(float2 x, float2 n){ return fract(x / n) * n;}
float3 emod3(float3 x, float3 n){ return fract(x / n) * n;}
/*** [BEGIN] SDF FUNCTIONS *******************************************/
float sdCircle(float2 p, float r) {
return length(p) - r;
}
float sdBox(float2 p, float2 b) {
float2 d = fabs(p) - b;
return length(fmax(d, (float2)(0.0f, 0.0f))) + fmin(fmax(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 = fmin(
fmin(
(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 = mix(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, fract(6.0f*q.y+sin(33.0f*q.x*q.y))), 3.0f)*(float3)(1.0f);
return pow(smoothstep(0.0f, 1.0f, fract(6.0f*q.y+sin(33.0f*q.x))), 3.0f)*(float3)(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(emod3(x3, (float3)(3.0f)), (float3)(1.0f)); // emod
float3 y = step(emod3(y3, (float3)(3.0f)), (float3)(1.0f)); // emod
float cell = step(emod(fragCoord.x*CELL_WIDTH+floor(fragCoord.y*CELL_WIDTH), 2.0f), 1.0f); // emod
return mix(x, y, cell);
}
float3 paintCheckerboard(float2 q, float h) // NEW
{
int2 cell = convert_int2(floor(q / h));
return (float3)((cell.x + cell.y) & 1);
}
float3 splatBlobRows(float2 p, float3 colBG, float3 colBlob, float HX, float HY, float dY, float blobRadius, int 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.5453f * ySnap) * 437.5453f; // 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 mix(colBG, colBlob, reg);
}
// 08088c,1f1170,1197f7,6119bf,ebfbfc,ffd4eb
__constant float3 paletteBluePurple[6] =// NEW
{
(float3)(0.0313725, 0.0313725, 0.5490196),
(float3)(0.1215686, 0.0666667, 0.4392157),
(float3)(0.0666667, 0.5921569, 0.9686275),
(float3)(0.3803922, 0.0980392, 0.7490196),
(float3)(0.9215686, 0.9843137, 0.9882353),
(float3)(1.0000000, 0.8313725, 0.9215686),
};
//------------------------------------------------------------------------------
// Palette 1: Coolors example (brown blue)
// 1a344d,1e8eb6,fc9843,5a1800,f7f5f6
//------------------------------------------------------------------------------
__constant 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
//------------------------------------------------------------------------------
__constant 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, __constant float3 *palette)// 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, 0);
col = splatBlobRows(p, col, palette[1], HX, HY, dY, 2.f * HX, 1);
col = splatBlobRows(p, col, palette[2], HX, HY, dY, HX / 2., 0);
col = splatBlobRows(p, col, palette[3], HX, HY, dY, HX / 3., 1);
col = splatBlobRows(p, col, palette[4], HX, HY, dY, HX / 6., 0);
col = splatBlobRows(p, col, palette[5], HX, HY, dY, HX / 7., 1);
return col;
}
float3 paintBlobs(float2 p) // NEW
{
return paintBlobsWithPalette(p, paletteGoldenMeadow);
}
/*****************************************************************************/
/*** [ END ] COLOR & TEST IMAGES ***********************/
/*****************************************************************************/
/*** [BEGIN] UTILITIES *******************************************/
__attribute__((overloadable)) 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 fract(sin(n)*43758.5453f); }
/*** [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) * fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(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) * fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(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 = fmax(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(fabs(Y) < 0.0001f) return M_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 = fabs(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 = fabs(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 = fabs(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 = fabs(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 = fabs(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 = fabs(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(fabs(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(fabs(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(fabs(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 && fabs(etaInf_p) < etaMax) {
// EXTRAPOLATE INWARD ON y=0 PARAMETERIZATION'S SEAM:
bool inWedge = (acos(fabs(normalize(p).x)) < 0.05f);
float xCrease = sign(dt) * (0.5f + fabs(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 *= fabs(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 = fabs(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); // 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*fabs(dt));
if(insertRemoveTine && applyPatch) {
if(fabs(p0.y) < cuspH && (sign(dt)*p0.x < 0.f || sign(dt)*p0.x > fabs(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 *= fabs(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);
// 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 * fabs(fract(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);
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);
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(fabs(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);
// 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 k=0; k<nPasses; k++) {
float yOffset = (float)(k)*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
float4 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(fmax(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 (float4)(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 = fmin(distLeft, 0.10f * fmax(eps, R)); // adaptive stepsize (spatially continuous errors)
distLeft -= dL;
// MIDPOINT STEP: dp = K(vmid) db
float2 db = dL * dir; // Δbrush
float4 K = getMixwellBrushMatrix(v, eps); // K(v)
float2 dp = (float2)(dot(K.xy, db), dot(K.zw, db)); // Euler approx, dp = K(v) db
float2 vmid = v + 0.5f * (dp - db); // v' = p' - b' = v + 0.5*(dp - db)
float4 Kmid = getMixwellBrushMatrix(vmid, eps); // K(vmid)
float2 dpmid = (float2)(dot(Kmid.xy, db), dot(Kmid.zw, 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(float2 p, float r, float2 a, 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] = {(float2)(p0), (float2)(p1), (float2)(p2), (float2)(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, fabs(sdf));
// JAGGED CASE (as rdrift superposition - symmetric but nonphysical)
float2 u = (float2)(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 and can produce a runaway loop that hangs the GPU. 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);
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 OpenCL for Houdini. Released under the MIT license.