Mixwell.glsl
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 ////////
///////////////////////////////////////////////////////////////////////
const float PI = 3.14159265358979323846; // π
const float PI_2 = 1.57079632679489661923; // π/2
// Constant that determines cutoff distance (as multiple of tine radius) beyond which rdSegment is zero.
const float RDSEGMENT_MASK_R_FACTOR = 10.;
/*** [BEGIN] SDF FUNCTIONS *********1**********************************/
float sdCircle( vec2 p, float r )
{
return length(p) - r;
}
float sdBox( in vec2 p, in vec2 b )
{
vec2 d = abs(p)-b;
return length(max(d,0.0)) + min(max(d.x,d.y),0.0);
}
// Unsigned distance from the infinite line specified by a ray origin (ro) and direction (unit rd).
float udLine(vec2 p, vec2 ro, vec2 rd) {
vec2 v = p-ro;
return length( v - dot(v,rd)*rd/dot(rd,rd) );
}
float sdSegment( in vec2 p, in vec2 a, in vec2 b )
{
vec2 pa = p-a, ba = b-a;
float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
return length( pa - ba*h );
}
float sdTriangle( in vec2 p, in vec2 p0, in vec2 p1, in vec2 p2 )
{
vec2 e0 = p1-p0, e1 = p2-p1, e2 = p0-p2;
vec2 v0 = p -p0, v1 = p -p1, v2 = p -p2;
vec2 pq0 = v0 - e0*clamp( dot(v0,e0)/dot(e0,e0), 0.0, 1.0 );
vec2 pq1 = v1 - e1*clamp( dot(v1,e1)/dot(e1,e1), 0.0, 1.0 );
vec2 pq2 = v2 - e2*clamp( dot(v2,e2)/dot(e2,e2), 0.0, 1.0 );
float s = sign( e0.x*e2.y - e0.y*e2.x );
vec2 d = min(min(vec2(dot(pq0,pq0), s*(v0.x*e0.y-v0.y*e0.x)),
vec2(dot(pq1,pq1), s*(v1.x*e1.y-v1.y*e1.x))),
vec2(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 :/
vec4 jet(float x) {
float r = clamp((x < 0.7) ? (4.0 * x - 1.5) : (-4.0 * x + 4.5), 0.0, 1.0);
float g = clamp((x < 0.5) ? (4.0 * x - 0.5) : (-4.0 * x + 3.5), 0.0, 1.0);
float b = clamp((x < 0.3) ? (4.0 * x + 0.5) : (-4.0 * x + 2.5), 0.0, 1.0);
return vec4(r, g, b, 1.0);
}
vec3 rgb(int r, int g, int b) {
return vec3(float(r), float(g), float(b))/255.;
}
vec3 quantizeColor(vec3 c, float n, bool bands) {
vec3 res = floor(c*n)/n;
if(bands) res *= smoothstep(-0.0001, +0.0001, length(c-res)-0.005);
return res;
}
// Paint texture
vec3 paintChocolate(vec2 q) {
float theta = 0.3;
float c = cos(theta);
float s = sin(theta);
vec2 p = vec2(c*q.x - s*q.y, s*q.x + c*q.y);
float k0 = 45.;
float k1 = 29.;
float phi0 = 3.;
float phi1 = 2.;
vec3 col0 = vec3(0.6862745098, 0.4 , 0.1490196078); //rgb(175, 102, 38);
vec3 col1 = vec3(0.6901960784, 0.9411764706, 0.9921568627); //rgb(176, 240, 253);
vec3 colChocolate = vec3(0.1490196078, 0.0862745098, 0.0901960784); //rgb(38,22,23);
float splat0 = smoothstep(0.5, 0.64, sin(k0*p.x*(1.-0.153*p.x)+phi0)*sin(k0*p.y*(1.+0.6*p.x)+phi0)*(1.-0.2*sin(12.*p.x*p.x)));
float splat1 = smoothstep(0.5, 0.64, sin(k1*p.x*(1.+0.230*p.x)+phi1)*cos(k1*p.y*(1.-0.1*p.x)+phi1)*(1.-0.2*sin(12.*p.x*p.y)));
vec3 col = mix(colChocolate + splat0*col0, col1, splat1);// + splat1*col1;
// TEST SQUARE:
//col = vec3(1.)* smoothstep(0.1, 0.11, q.x) * smoothstep(0.2, 0.19, q.x) * smoothstep(0.1, 0.11, q.y) * smoothstep(0.2, 0.19, q.y);
return col;
}
vec3 paintStripes(vec2 q) {
//return pow(smoothstep(0., 1., fract(6.*q.y+sin(33.*q.x*q.y))),3.)* vec3(1);// < 0.5 ? vec3(1) : vec3(0);
return pow(smoothstep(0., 1., fract(6.*q.y+sin(33.*q.x))),3.)* vec3(1);// < 0.5 ? vec3(1) : vec3(0);
}
// Adapted from https://webgl-operate.org/examples/canvassize-example.html
vec3 paintTestImage(vec2 fragCoord) {
const float CELL_WIDTH = 1.0 / 32.0;
vec3 x3 = vec3(fragCoord.x) + vec3(0.0, 1.0, 2.0);
vec3 y3 = vec3(fragCoord.y) + vec3(0.0, 1.0, 2.0);
vec3 x = step(mod(x3, vec3(3.0)), vec3(1.0));
vec3 y = step(mod(y3, vec3(3.0)), vec3(1.0));
float cell = step(mod(fragCoord.x*CELL_WIDTH + floor(fragCoord.y*CELL_WIDTH), 2.0), 1.0);
return mix(x, y, cell);
}
vec3 paintCheckerboard(vec2 q, float h) {
return vec3(mod( dot(mod ( floor(q/h), 2.),vec2(1)), 2.));
}
vec3 paintTestSquare(vec2 q) { // NEW
return vec3(1.0) * smoothstep(0.1, 0.11, q.x) * smoothstep(0.2, 0.19, q.x) *
smoothstep(0.1, 0.11, q.y) * smoothstep(0.2, 0.19, q.y);
}
vec3 splatBlobRows(vec2 p, vec3 colBG, vec3 colBlob, float HX, float HY, float dY, float blobRadius, bool isOddRow) { // NEW
float Yoffset = isOddRow ? 0.5*HY : 0.;
// 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. + ySnap; // nat freq for variation
// Mission: Find our blob i:
float x0 = HX*floor(p.x/HX);
float x1 = x0 + HX; // ceil
vec2 p0 = vec2(x0, dY*sin(F*x0));
vec2 p1 = vec2(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);
vec2 pC = vec2(xC, yC);
vec2 pL = vec2(xL, yL);
vec2 pR = vec2(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./DC - 1./DL - 1./DR;
// Implicit blob size hint:
f -= 1./(blobRadius*blobRadius);
float reg = smoothstep(HX, 2.*HX, max(0., f));
return mix(colBG, colBlob, reg);
}
const vec3 paletteBluePurple[6] = vec3[6]( // NEW
vec3( 8., 8., 140.)/255.0,
vec3( 31., 17., 112.)/255.0,
vec3( 17., 151., 247.)/255.0,
vec3( 97., 25., 191.)/255.0,
vec3(235., 251., 252.)/255.0,
vec3(255., 212., 235.)/255.0
);
//------------------------------------------------------------------------------
// Palette 1: Coolors example (brown blue)
// 1a344d,1e8eb6,fc9843,5a1800,f7f5f6
//------------------------------------------------------------------------------
const vec3 paletteCoolors1[6] = vec3[6]( // NEW
vec3(0.1019608, 0.2039216, 0.3019608),
vec3(0.3529412, 0.09411765, 0.0),
vec3(0.1176471, 0.5568628, 0.7137255),
vec3(0.9882353, 0.5960785, 0.2627451),
vec3(0.9686275, 0.9607843, 0.9647059),
vec3(0.9686275, 0.9607843, 0.9647059)
);
//------------------------------------------------------------------------------
// Palette: fb6107,f3de2c,7cb518,5c8001,fbb02d
// Names: Blaze Orange, Golden Glow, Lime Moss, Forest Moss, Sunflower Gold
//------------------------------------------------------------------------------
const vec3 paletteGoldenMeadow[6] = vec3[6]( // NEW
vec3(0.9843137, 0.3803922, 0.0274510), // fb6107
vec3(0.3607843, 0.5019608, 0.0039216), // 5c8001
vec3(0.9843137, 0.6901961, 0.1764706), // fbb02d
vec3(0.4862745, 0.7098039, 0.0941176), // 7cb518
vec3(0.9529412, 0.8705882, 0.1725490), // f3de2c
vec3(0.9529412, 0.8705882, 0.1725490) // f3de2c
);
vec3 paintBlobsWithPalette(vec2 p, vec3 palette[6]) { // NEW
float HX = 0.25; // X CELL SIZE
float HY = 1.; // ROW SPACING
float dY = HX/10.; // Y POSITION VARIATION
float blobRadius = 2.*HX; // BLOB SIZE HINT FOR SMALLER BLOBS
vec3 col = vec3(0.); // Background
col = splatBlobRows(p, col, palette[0], HX, HY, dY, 2.*HX, false);
col = splatBlobRows(p, col, palette[1], HX, HY, dY, 2.*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;
}
vec3 paintBlobs(vec2 p) { // NEW
return paintBlobsWithPalette(p, paletteGoldenMeadow);
}
/*****************************************************************************/
/*** [ END ] COLOR & TEST IMAGES ***********************/
/*****************************************************************************/
/*** [BEGIN] UTILITIES *******************************************/
vec2 rotate(vec2 v, float theta) {
float c = cos(theta);
float s = sin(theta);
return vec2(
c * v.x - s * v.y,
s * v.x + c * v.y
);
}
// Returns (-y,x) given v=(x,y).
vec2 rot90(vec2 v) {
return vec2(-v.y, v.x);
}
// Normalizes p if ‖p‖<1.
vec2 projectOutsideUnitDisk(vec2 p) {
float R2 = p.x*p.x + p.y*p.y;
return (R2 >= 1.) ? p : p/sqrt(R2);
}
// Simple pseudorandom float-to-float[0,1] hash.
float iqhash(float n) { return fract(sin(n)*43758.5453); }
/*** [END] UTILITIES *******************************************/
/*****************************************************************************/
/*** [BEGIN] CARLSON ITERATIONS AND ELLIPTIC FUNCTIONS ***********************/
/*****************************************************************************/
const int CARLSON_MAX_ITER = 15; // Maximum number of iterations
const float CARLSON_TOL = 0.001; // Tolerance for convergence
// Carlson's RF function [Carlson 1995, p.16] (as .x) and the final "n" iteration# (as .y).
vec2 CarlsonRFn(float x, float y, float z) {
float A0 = (x + y + z) * 0.3333333333333333;;
float An = A0;
float Q = pow(3.0 * CARLSON_TOL, -1.0/6.0) * max(max(abs(A0-x), abs(A0-y)), abs(A0-z));
float xn = x;
float yn = y;
float zn = z;
float n4 = 1.0;
int n = 0;
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.25;
xn = (xn + lambda) * 0.25;
yn = (yn + lambda) * 0.25;
zn = (zn + lambda) * 0.25;
Q *= 0.25;
n4 *= 4.0; // 4^n for break or CARLSON_MAX_ITER
}
float X = (A0 - x)/(n4 * An);
float Y = (A0 - y)/(n4 * An);
float Z = -X - Y;
float E2 = X*Y - Z*Z; // X*Y + Y*Z + Z*X;
float E3 = X*Y*Z;
float RF = (1. - (E2/10.) + (E3/14.) + (E2*E2/24.) - (3./44.*E2*E3)) / sqrt(An);
return vec2(RF,n);
}
// Carlson's RD function [Carlson 1995, p.20] (as .x) and the final "n" iteration# (as .y).
vec2 CarlsonRDn(float x, float y, float z) {
float A0 = (x + y + 3.0*z)/5.0;
float An = A0;
float Q = pow(CARLSON_TOL/4.0, -1.0/6.0) * max(max(abs(A0-x), abs(A0-y)), abs(A0-z));
float xn = x;
float yn = y;
float zn = z;
float Sum = 0.0;
float n4 = 1.0;
int n = 0;
for (n=0; n<CARLSON_MAX_ITER; n++) {
float lambda = sqrt(xn * yn) + sqrt(xn * zn) + sqrt(yn * zn);
Sum += 1.0 / ( n4 * sqrt(zn) * (zn + lambda) );
if (Q <= CARLSON_TOL * An) {
break;
}
An = (An + lambda) * 0.25;
xn = (xn + lambda) * 0.25;
yn = (yn + lambda) * 0.25;
zn = (zn + lambda) * 0.25;
Q *= 0.25;
n4 *= 4.0; // 4^n for break or CARLSON_MAX_ITER
}
float X = (A0 - x)/(n4 * An);
float Y = (A0 - y)/(n4 * An);
float Z = -(X+Y)/3.;
float XY = X*Y;
float Z2 = Z*Z;
float E2 = XY - 6.*Z2;
float E3 = (3.*XY - 8.*Z2)*Z;
float E4 = 3.*(XY - Z2)*Z2;
float E5 = XY*Z2*Z;
float Expansion = (1. - 3./14.*E2 + E3/6. + 9./88.*E2*E2 - 3./22.*E4 - 9./52.*E2*E3 + 3./26.*E5);
float RD = Expansion/(n4 * An*sqrt(An)) + 3.*Sum;
return vec2(RD,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.0 - m*sinPhi*sinPhi, 1.0); /// 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.0 - m*sin2Phi, 1.0)
- (m*sin3Phi/3.0) * CarlsonRD(cos2Phi, 1.0 - m*sin2Phi, 1.0);
}
// Mathematica-style Incomplete Elliptic Integrals of both the First and Second Kinds,
// EllipticF(φ,m) and EllipticE(φ,m), respectively, with m=k² parameter.
// Computed in terms of Carlson's RF & RD functions [Carlson 1995, (4.5,4.6)].
// Faster than independent calls since shares a common CarlsonRF evaluation.
// Also returns the max number of Carlson RD | RF iterations, n.
// @return vec3(F,E,n)
vec3 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.0 - m*sinPhi2;
vec2 RFn = CarlsonRFn(x, y, 1.0);// (RF,nRF)
vec2 RDn = CarlsonRDn(x, y, 1.0);// (RD,nRD)
float n = max(RFn.y, RDn.y);// max(nRF,nRD)
float F = sinPhi * RFn.x;
float E = F - (m*sinPhi3/3.0) * RDn.x;
return vec3(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 vec2(F,E)
vec2 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)
vec2 mxwGradEtaInf(vec2 p) {
float X = p.x;
float Y = p.y;
float R2 = dot(p,p);
float R4 = R2*R2;
return vec2( 2.*X*Y/R4, 1. + 2.*Y*Y/R4 - 1./R2 );
}
// Signed η_∞ for point p=(X,Y) outside the unit circle.
float mxwEtaInf(vec2 p) {
float X = p.x;
float Y = p.y;
return Y - Y/(X*X + Y*Y);
}
// Signed η_∞ for point p=(X,Y) outside the unit circle.
float mxwEtaInf(float X, float 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(vec2 p) {
float X=p.x; float Y=p.y;
//if(abs(Y) < 0.0001) 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.0) ? ((sigma>=0.0) ? (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(vec2 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); // σ∈[-π/2,π/2]
sigma = mxwClampSigmaPosY(sigma); // Safety first.
float etaInf = Y - Y/R2;
float etaInf2 = etaInf*etaInf;
float etaSq4 = etaInf2 + 4.0;
float m = 4.0/etaSq4;
vec2 FE = EllipticFE(sigma,m);
float F = FE.x;
float E = FE.y;
float x = ( etaSq4*E - (etaSq4-2.0)*F )/(2.0*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.0;
float m = 4.0/etaSq4;
vec2 FE = EllipticFE(sigma,m);
float F = FE.x;
float E = FE.y;
float x = ( etaSq4*E - (etaSq4-2.0)*F )/(2.0*sqrt(etaSq4));
float cosS = cos(sigma);
float Y = 0.5*(eta + sqrt(eta2 + 4.0*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.
vec2 mxwPositionAtSigmaEta(float sigma, float eta) {
float etaAbs = abs(eta);
float eta2 = eta*eta;
float cosS = cos(sigma);
float fug = (eta2 + 4.0*cosS*cosS);// abs ∵ floating-point demons (!)
float Y = 0.5*(etaAbs + sqrt(eta2 + 4.0*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 vec2(X,Y);
}
// Maxwell arctan monotone preconditioner, M, and its derivative, D, at (t,η∞).
// Returns vec2(M,D).
vec2 mxwMD_arctan(float t, float eta) {//eta>0
// REUSE S(eta):
float eta2 = eta*eta;
float S = 0.5*( eta + (eta2 + 2.0)/sqrt(eta2 + 4.0) ); // Use abs(eta) to make it symmetric top/bottom
float twoOverPI = 0.6366197723675813430;
float M = +twoOverPI * atan(-t/S);
float D = -twoOverPI / (S + t*t/S);
return vec2(M, D);
}
// Maxwell tanha monotone preconditioner, M, and its derivative, D, at (t,η∞).
// Returns vec2(M,D).
vec2 mxwMD_tanha(float t, float eta) {//eta>0
// REUSE S(eta):
float eta2 = eta*eta;
float S = 0.5*( eta + (eta2 + 2.0)/sqrt(eta2 + 4.0) ); // Use abs(eta) to make it symmetric top/bottom
// Optimized magic o((⊙﹏⊙))o numbers
float a = 0.956764;
float s = 1.05782;
float c = 1.0 / (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 vec2(M, D);
}
// Maxwell monotone preconditioner's magic α* value (for η≥0).
float mxwMD_alphaStar(float eta) {
float n = eta;// ouch, lazy
float n2 = eta*eta;
float squirt = sqrt(n2+4.0);
float alpha = 2.0*n*squirt / (n2 + n*squirt + 2.0);
return alpha;
}
// Maxwell monotone preconditioner M, and its derivative, D, at specified (t,η∞).
// Returns vec2(M,D).
vec2 mxwMD(float t, float eta) {
float absEta = abs(eta);
vec2 arctanMD = mxwMD_arctan(t, absEta);
vec2 tanhaMD = mxwMD_tanha (t, absEta);
float alpha = mxwMD_alphaStar(absEta);
vec2 MD = alpha*arctanMD + (1.0-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; // ooh, lazy.
float n2 = eta*eta;
float cosS = cos(sigma);
float cos2 = cosS*cosS;
float squirt = sqrt(n2 + 4.0*cos2);
return -0.5 * (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:
vec2 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++) {// MD_arctan MD_axwell
float t = mxwTimeAtSigmaEta(sigma, etaAbs);
vec2 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;
sigma = mxwClampSigmaPosY(sigma); // Safety first
// ADAPTIVE NEWTON TOL:
float SIGMA_NEWTON_TOL = 0.0000002; // Newton solver tol for σ(t)
float tolExp = 5.0*smoothstep(0.01, 0.5, etaAbs);
SIGMA_NEWTON_TOL *= pow(2., tolExp);
if(etaAbs > 1.) SIGMA_NEWTON_TOL = 0.00001;
if(abs(dSigma) < PI_2 * SIGMA_NEWTON_TOL) break;
}
// REFLECT IF η<0: ensure π/2<|σ|<π
if(eta<0.0) {
if(sigma>=0.0)
sigma = PI - sigma;
else
sigma = -(PI + sigma);
}
return sigma;
}
// RAW MONOTONE-PRECONDITIONED NEWTON SOLVE FOR POSITION AT (t,η).
// POOR NEAR η≅0.
vec2 mxwPositionSolveAtTimeEta(float t, float eta) {
eta = (eta!=0.0) ? eta : 0.0001; // Please stay off the line folks.
if(abs(eta)<0.0001) eta = 0.0001 * sign(eta);
float sigma = mxwSigmaSolver(t, eta); // RETURNS σ∈(π,-π). Bad for large |t|.
vec2 p = mxwPositionAtSigmaEta(sigma, eta); // POOR NEAR η≅0.
return p;
}
// RAW MONOTONE-PRECONDITIONED NEWTON SOLVER FOR POSITION ADVECTED FOR DT TIME.
vec2 mxwAdvectPointRaw(vec2 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 eta~0.
}
// Advects point through unit-cylinder-flow but transitions to "pure translation" outside a
// maximum T/Dist range to avoid distortion near x-axis.
// Uses Maxwell time monotone-preconditioned Newton advection solver (MTMPNAS).
//
// 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.
vec2 mxwAdvectPointTCutoff(vec2 p, float dt, float tCutoff) {
float eta = mxwEtaInf(p);
float t = mxwTimeAt(p);// BUG: Huge values near eta~0
float tNew = t + dt;
vec2 pNew;
if(abs(tNew)<tCutoff) {
pNew = mxwPositionSolveAtTimeEta(tNew, eta); // BUG: HUGE VALUES OF tNew near η∞~0.
}
else {// Break advection into two steps: (1) near-field short-time σ(t) solve, and (2) far-field translation.
float tSafe = sign(tNew) * tCutoff;
float tSlide = tNew - tSafe;
pNew = mxwPositionSolveAtTimeEta(tSafe, eta);
pNew.x -= tSlide;
}
return pNew;
}
// Advects p for time dt through flow around unit circle.
vec2 mxwAdvectPointUnpatched(vec2 p, float dt) {
return mxwAdvectPointTCutoff(p, dt, 20.);// Hard-coded TCutoff value
}
vec2 mxwAdvectPointUnitPatched(vec2 p, float dt, bool applyPatch) {
float etaInf_p = mxwEtaInf(p); /// N.B.: Computer AFTER tine insertion.
// ADVECT TEXTURE AROUND UNIT CYLINDER:
vec2 pNew = vec2(p);
float etaMax = 0.05; // patch size
if(applyPatch && abs(etaInf_p) < etaMax)
{
//vec2 a = vec2(0.);
//vec2 b = vec2(dt, 0.);
//pNew = p - b + rdSegment(p, b, a, 1.0);
// EXTRAPOLATE INWARD ON y=0 PARAMETERIZATION'S SEAM:
bool inWedge = (acos(abs(normalize(p).x)) < 0.05);
float xCrease = sign(dt) * (0.5 + abs(dt));
bool inPatchRegion = (inWedge &&
( (dt<0. && (p.x>= 1.0 || p.x<xCrease)) ||
(dt>0. && (p.x<=-1.0 || p.x>xCrease)) ) );
if(inPatchRegion) {
// extrapolate from off-line y∈[-H,H]
float H = etaMax;
float dy = (p.y<0. ? -H : H);
vec2 pPerturb = vec2(p.x, dy );
vec2 pPerturbNew = mxwAdvectPointUnpatched(pPerturb, dt); // [1st advection eval]
vec2 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):
// float uy0 = u.y; float uy1 = u2.y;
// float a = (8.*uy0 - uy1)/(6.*H);
// float b = (uy1 - 2.*uy0)/(6.*H*H*H);
// float y = abs(p.y);
// float uy = a*y + b *y*y*y;
// u.y = uy;
//}
// MINOR CURVATURE CORRECTION (noticeable around unit-cyl flow, but not cusp brush)
// Avoid for cusp brush since uses an extra sample.
if(false) {// FANCY u.x CURVATURE CORRECTION:
vec2 pPertur2 = vec2(p.x, dy*2.);
vec2 pPertur2New = mxwAdvectPointUnpatched(pPertur2, dt); // [2nd advection eval in patch]
vec2 u2 = (pPertur2New - pPertur2); // for u.x curvature correction
// PARABOLIC MODEL FOR u.x(y) CORRECTION: (call ux -> x)
// x(y) = a y² + b
// x1=x(h) & x2=x(2h)
// → a=(x2-x1)/(3h²) & b=(4 x1 - x2)/3
{
float ux1 = u.x;
float ux2 = u2.x;
float a = (ux2-ux1)/(3.*H*H);
float b = (4.*ux1 - ux2)/3.;
float y = p.y;
u.x += (a*y*y + b - u.x);
}
}
pNew = p + u;
return pNew;
}// PATCH END
}// ELSE NOT IN PATCH REGION:
// MONOTONE-PRECONDITIONED SOLVER:
pNew = mxwAdvectPointUnpatched(p, dt); // [single advection eval]
return pNew;
}
// Expands p to accommodate 2D incompressible insertion of a tine of radius r at the origin.
vec2 mxwApplyTineInsertionMap(vec2 p, float r) {
if(dot(p,p)==0.) p = 0.000001*vec2(1.); // avoids R=0 and 0/0 result.
float R = length(p); // current radius
float S = sqrt(R*R + r*r);
return S/R*p;
}
// Contracts p to accommodate 2D incompressible removal of a tine of radius r at the origin.
vec2 mxwApplyTineRemovalMap(vec2 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;
}
// Maxwell unit-cylinder advection solver with tine insertion/removal and singularity patch.
vec2 mxwAdvectPointUnit(vec2 p, float dt, bool insertRemoveTine, bool applyPatch)
{
vec2 p0 = p;// initial position
// SETUP CUSP PATCH:
bool cuspPatch = false;
float cuspH = min(0.02, 0.005*abs(dt));
if(insertRemoveTine && applyPatch) {// PATCH THIS CUSP BRUSH:
//cuspH = min(cuspH, abs(p0.x));//wedge
if(abs(p0.y)<cuspH && (sign(dt)*p0.x<0.0 || sign(dt)*p0.x > abs(dt)*0.8) ) {
float dy = (p0.y<0. ? -cuspH : cuspH);
//dy = 0.995*dy + 0.005*p0.y;// breaks-up constant-extrap. color banding
p.y = dy; // perturb compute location
cuspPatch = true; // apply blend at end to u.y
}
}
vec2 pPreCusp = p;// for computing displacement, u
// INSERTION: INFLATE COORDS for R=1+epsRegularization circle:
float epsR = (applyPatch ? 0.0005 : 0.0); // pull out a little more paint when rd-removing tine
p = (insertRemoveTine ? mxwApplyTineInsertionMap(p, 1.+epsR) : p); // rd end pos
// ADVECTION AROUND PATCHED UNIT CYLINDER:
vec2 pNew = mxwAdvectPointUnitPatched(p, dt, applyPatch);
// IF REMOVE TINE, DEFLATE COORDS:
pNew = (insertRemoveTine ? mxwApplyTineRemovalMap(pNew, 1.0) : pNew); // rd start pos
// PATCH CUSP:
if(cuspPatch) {
vec2 u = pNew - pPreCusp; //p0;
u.y *= abs(p0.y)/cuspH;
pNew = p0 + u; // filtered u.y
}
return pNew;
}
// Advect point p around a radius-R cylinder at the origin for flow distance L in -ve X direction.
vec2 mxwAdvectPointScaled(vec2 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:
vec2 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) {
float y = height/r;
float dx = 0.; //drift = 5.*r/(0.2 + pow(eta,3.));
float eps = 0.002;
y = sqrt(y*y + eps*eps); // y_eps
float y2 = y *y;
float y4 = y2*y2;
float y6 = y4*y2;
float y8 = y6*y2;
float y11 = y8*y2*y;
if(y > 6.50416858646776)
{
dx = -((PI*(6615. - 1960.*y2 + 600.*y4 - 192.*y6 + 64.*y8))/(256.*y11));
}
else if (y > 0.8220844420096408)
{
float dy = y - 2.3;
float dy2 = dy*dy;
float dy4 = dy2*dy2;
dx = (-0.0413839 - 0.0937878*dy - 0.043355*dy2 - 0.0068946*dy2*dy + (1.77539e-6)*dy4 - (1.01973e-7)*dy4*dy) / (1.00000000000000 + 3.26035*(dy) + 3.66495*dy2 + 2.08615*dy2*dy + 0.665761*dy4 + 0.115772*dy4*dy + 0.00875686*dy4*dy2);
}
else
{
dx = -0.039720771 - 0.16369764*y2 + 0.0076619254*y4 - 0.00096477622*y6 + 0.00015811568*y8
+ (0.50000000 + 0.093750000*y2 - 0.0073242188*y4 + 0.0010681152*y6 - 0.00018775463*y8)*log(y);
}
return 2.0 * r * dx; // UNDO UNIT-RADIUS XFORM & *2 for full drift.
}
// Returns rdLine1DS series implementation of 1D drift. Provided for "backwards convenience." // NEW
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
vec2 rdSegmentM(vec2 p, float r, vec2 a, vec2 b) {
if(length(a-b)==0.) return vec2(0.);
// Reverse of a→b motion is a←b:
// So shift b to the origin (start):
p -= b;
vec2 M = a - b; // Reverse motion vector
// Rotate to x-axis flow:
float L = length(M); // assume dt=L/r>0 always
vec2 newXDir = normalize(M);
//vec2 newYDir = vec2(-newXDir.y, newXDir.x);
float cs = newXDir.x;
float sn = newXDir.y;
mat2 Q = mat2(cs, -sn, +sn, cs); //rotation of X to B
mat2 QT = mat2(cs, +sn, -sn, cs); //rotation of B to X
if(dot(p,p)==0.) return M;// point at cusp tip.
vec2 p0 = p;
p = Q*p; // Rotate to x-axis flow orientation
bool insertRemoveTine = true;
bool applyPatch = true;
vec2 pNew = mxwAdvectPointScaled(p, L, r, insertRemoveTine, applyPatch);
pNew = QT*pNew; // Unrotate back to B flow orientation
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.002;
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.0/y2;
vec2 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.)*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: vec2 p - 2D point.
// Output: vec2 v - Velocity vector in World frame.
vec2 cylFlowVelWorld(vec2 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.0 / (R2 * R2); // Inverse R⁴
return vec2((X2 - Y2) * invR4, 2.0 * 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: vec2 p - 2D point. Output: vec2 - Velocity vector in Body frame.
vec2 cylFlowVelBody(vec2 p) {
return cylFlowVelWorld(p) - vec2(1.,0.);
}
// 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: vec2 p - 2D point. Unit velocity of cylinder, vc.
// Output: vec2 - Velocity vector in World frame.
// NOTE: Matches cylFlowVelWorld(p) when vc=(1,0).
//vec2 cylFlowVelWorldArbDir(vec2 p, vec2 vc) {
// vec2 d = normalize(vc); // unit direction of cylinder motion (SLOTH)
// mat2 R = mat2(d.x, -d.y,
// d.y, d.x); // Rotates (1,0) to d
// mat2 RT = transpose(R); // Rotates d to (1,0)
//
// // PULL BACK TO REFERENCE FRAME (flow in -X direction)
// vec2 pRef = R*p;
// vec2 vRef = cylFlowVelWorld(pRef);
// vec2 v = RT*vRef; // PUSH FWD TO WORLD ORIENT
// return v;
//}
vec2 cylFlowVelWorldArbDir(vec2 p, vec2 vc) {
vec2 d = normalize(vc); // unit direction of cylinder motion (SLOTH)
//mat2 R = mat2(d.x, -d.y, d.y, d.x); // Rotates (1,0) to d
//mat2 RT = transpose(R); // Rotates d to (1,0)
// PULL BACK TO REFERENCE FRAME (flow in -X direction)
vec2 pRef = vec2(d.x*p.x + d.y*p.y, -d.y*p.x + d.x*p.y); // R*p
vec2 vRef = cylFlowVelWorld(pRef);
// PUSH FWD TO WORLD ORIENT
return vec2(d.x*vRef.x - d.y*vRef.y, d.y*vRef.x + d.x*vRef.y); // RT*vRef
}
// 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.; // big enough to avoid negative drift
vec2 p0 = vec2(Xmax, y); // Particle init at (0,y) --> (pureTranslation-dx, ~y)
vec2 p = p0;
int i = 0;
float time = 0.0;
while(p.x>-Xmax && i<3000) {
float dt = 0.01 * length(p); // adaptive stepsize
vec2 vp = cylFlowVelBody(p); // initial velocity; (-1,0) at ∞
vec2 pmid = p + vp*dt*0.5; // midpoint position
vec2 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.; // big enough to avoid negative drift
vec2 p0 = vec2(0., y);// Particle init at (0,y) --> (-dx, ~y)
vec2 c = vec2(-Xmax,0.);// Cylinder position
vec2 vc = vec2(+1.0, 0.);// Cylinder velocity (constant for 1D drift calc)
vec2 p = p0;
int i = 0;
float time = 0.0;
while(c.x<Xmax && i<3000) {
vec2 q = p - c; // relative particle position
float dt = 0.01 * length(q); // adaptive stepsize
vec2 vp = cylFlowVelWorld(q); // initial velocity; (-1,0) at ∞
vec2 pmid = p + vp*dt*0.5; // midpoint position
vec2 cmid = c + vc*dt*0.5; // midpoint cylinder position
vec2 qmid = pmid-cmid; // midpoint relative particle position
vec2 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.
vec2 rdLine(vec2 p, float r, vec2 dir)
{
dir = normalize(dir);
vec2 n = vec2(-dir.y, +dir.x);
float height = dot(p,n);
vec2 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(fract(y/hy-0.5)-0.5);
}
// Nonpareil pattern.
vec2 rdNonpareil(vec2 p, float combR, vec2 combDir, float combGap, int nPasses) {
vec2 dir = normalize(combDir);
vec2 n = rot90(dir);
float y = dot(p,n);
float passGap = float(nPasses)*combGap; // wider spacing between passes for nicer falloff + overlap
vec2 rdf = vec2(0.);
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;
}
// Nonpareil RDF 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).
vec2 rdNonpareilNoisy(vec2 p, float r, vec2 dir, float combGap, float combGapNoise) {
dir = normalize(dir);
vec2 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.
vec2 rdf = vec2(0.);
for(int k=k0-dkMax; k<=k0+dkMax; k++) {
float noisek = 2.*iqhash(float(k))-1.;// ∈[-1,1)
float yk = combGap * float(k) + combGapNoise*noisek;
////if(iqhash(noisek*PI) > 0.12) // Optional tine dropout (skip for now)
rdf += rdLine1DS(abs(y-yk),r) * dir;
}
return rdf;
}
// Gel-Git pattern.
vec2 rdGelGit(vec2 p, float combR, vec2 combDir, float combGap, int nPasses) {
vec2 dir = normalize(combDir);
vec2 n = rot90(dir);
float y = dot(p,n);
float passGap = float(nPasses)*combGap*2.0; // wider spacing between passes for nicer falloff + overlap
vec2 rdf = vec2(0.);
// GIT ("left")
for(int k=0; k<nPasses; k++) {
float yOffset = (0.5 + float(k))*combGap*2.0;
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.0;
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, vec3 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).
vec2 waveUnitTangent(float X, vec3 W) {
float Df = W.x * W.y * cos(W.y * X + W.z); // f'(X)
vec2 T = vec2(1.0, 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: vec3 W=(A,K,P) describes the partial wave A*sin(K*X+P)
// rightComb: true passes right (+x), false passes left (-x) direction.
vec2 rdSineT(vec2 p, float r, float waveY0, in vec3 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.0; // now X offset
float y0 = p.y; // INIT POS: p=(0,y0)
// WAVE PARAMS:
W = vec3(waveAmp, waveK, wavePhase); // scaled wave params
float rdir = (rightComb ? -1.0 : +1.0);// reverse x-dir for combing.
// ADVECT PARTICLE at p AS CYLINDER MOVES ON SINEWAVE FROM (x=-Xmax to +Xmax).
float Xmax = 50.; // big enough to avoid negative drift
vec2 p0 = vec2(0., y0); // Particle init at (0,y0) --> p
vec2 c = vec2(-Xmax*rdir, 0.);// Cylinder position (bogus y value)
vec2 vc; // Cylinder velocity (unit vector)
int i = 0;
float time = 0.0;
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)
vec2 q = p - c; // relative particle position
float dt = 0.05; // HQ:0.01 SQ:0.05
vec2 vp = cylFlowVelWorldArbDir(q,vc); // particle velocity
vec2 pmid = p + vp*dt*0.5; // midpoint position
vec2 cmid = c + vc*dt*0.5; // midpoint cylinder position
vec2 vcmid = waveUnitTangent(cmid.x,W) *rdir;
vec2 qmid = pmid - cmid; // midpoint relative particle position
vec2 vmid = cylFlowVelWorldArbDir(qmid,vcmid); // midpoint velocity
// FWD EULER:
//p += vp * dt; // position update
//c += vc * dt; // advance cylinder
// MIDPOINT METHOD:
p += vmid * dt;
c += vcmid * dt;
p = projectOutsideUnitDisk(p-c) + c; // Avoid losing particles inside disk
time += dt;
i++;
}
vec2 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) {
vec2 rdrift = rdSineT(vec2(0.,height), r, 0.0, vec3(0.0, 1.0, 0.0), true);
return rdrift.x;
}
// Combs a→b (performs insert/remove flows at begin/end)
// Time-stepped implementation (T).
vec2 rdSegmentT(vec2 p, float r, vec2 a, vec2 b) {
// UNIT-RADIUS TRANSFORM:
float invR = 1.0/r;
p *= invR;
a *= invR;
b *= invR;
// REVERSE COMB SETUP: b→a
vec2 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.
vec2 p0 = p; // Particle init (USED TO GET RevDrift: p-p0)
vec2 c = b; // Cylinder init
vec2 vc = normalize(rdir); // Cylinder unit velocity
int i = 0;
float time = 0.0;
// IF REMOVE TINE, INFLATE COORDS:
//p = applyTineInsertionMap(p, R ); // BODY FRAME
p = b + mxwApplyTineInsertionMap(p-b, 1.); // WORLD FRAME
// Reverse-advect p as cylinder moves from b→a:
float tLeft = L;
while(tLeft > 0.000) {
vec2 q = p - c; // relative particle position
float dt = 0.2;//max(0.01, 0.10*min(1.0, length(p-b)));
dt *= max(0.05, dist2Segment); // adaptive stepsize (near segment)
if(dt > tLeft) dt = tLeft;
// TODO: should depend on L, r.
vec2 vp = cylFlowVelWorldArbDir(q,vc); // particle velocity
vec2 pmid = p + vp*dt*0.5; // midpoint position
vec2 cmid = c + vc*dt*0.5; // midpoint cylinder position
vec2 qmid = pmid - cmid; // midpoint relative particle position
vec2 vmid = cylFlowVelWorldArbDir(qmid,vc); // midpoint velocity
// FWD EULER:
//p += vp * dt; // position update
//c += vc * dt; // advance cylinder
// 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.); // WORLD FRAME; unit tine
vec2 rdrift = p-p0; // reverse drift
return rdrift * r; // UNDO UNIT-RADIUS TRANSFORM
}
// Mixwell brush matrix for relative position, v, and radius eps. // NEW
mat2 getMixwellBrushMatrix(vec2 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 = inversesqrt(s);
float invS = invSqrtS*invSqrtS;
float invD = invS*invSqrtS;
float invR = inversesqrt(max(R2, 1e-20)); // 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.0 - R*(R2 + 2.0*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 mat2(Afd + xx, xy, xy, Afd + yy);
}
// Reverse drift through Mixwell flow of radius-eps brush moving a→b. Computed using adaptive midpoint integration. // NEW
vec2 rdMixwellBrushAdaptiveMidpoint(vec2 p0, float eps, vec2 a, vec2 b) {
// Step brush position backwards from b to a, advecting p to compute rev-drift (p-p0):
vec2 uBrush = a - b; // rev brush step
float L = length(uBrush);
if (L <= 1e-20) return vec2(0.);
vec2 dir = uBrush/L; // rev brush dir
float distLeft = L; // brush motion left
vec2 p = p0; // reverse-drifted point position (init: p0)
while (distLeft > 0.) {
vec2 v = p - b; // rel pos from end brush goal
float R = length(v);
float dL = min(distLeft, 0.10*max(eps, R)); // adaptive stepsize (spatially continuous errors)
distLeft -= dL;
// MIDPOINT STEP: dp = K(vmid) db
vec2 db = dL*dir; // Δbrush
mat2 K = getMixwellBrushMatrix(v, eps); // K(v)
vec2 dp = K*db; // Euler approx, dp = K(v) db
vec2 vmid = v + 0.5*(dp - db); // v' = p' - b' = v + 0.5*(dp - db)
mat2 Kmid = getMixwellBrushMatrix(vmid, eps); // K(vmid)
vec2 dpmid = 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
vec2 rdSegmentMBrush(vec2 p, float r, vec2 a, vec2 b) {
return rdMixwellBrushAdaptiveMidpoint(p, r, a, b);
}
// Combs a→b (performs insert/remove flows at begin/end).
vec2 rdSegment(vec2 p, float r, vec2 a, vec2 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 vec2(0.);
float mask = 1. - smoothstep(0.5*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
//if(d < 0.005f*r) // OPTIONAL: Timestep on-line singularities
// return mask * rdSegmentT(p, r, a, b); // Time-stepped advection solver
//else
// return mask * rdSegmentM(p, r, a, b); // Maxwell advection solver
}
// Combs a→b (performs insert/remove flows at begin/end), with reverse option.
// @param reverse - Swaps a↔b.
vec2 rdSegmentRev(vec2 p, float r, vec2 a, vec2 b, bool reverse) {
vec2 head = (reverse ? b : a);
vec2 tail = (reverse ? a : b);
return rdSegment(p, r, head, tail);
}
vec2 rdTriangle(vec2 p, float r, vec2 p0, vec2 p1, vec2 p2)
{
vec2 pInit = p; // copy to compute rdrift
vec2 P[4] = vec2[4] ( p0, p1, p2, p0 );
float sdf = sdTriangle(p, p0, p1, p2);
float sdfMax = 5.*r;
if(sdf > sdfMax) return p-pInit;
// Masking factor to avoid discontinuity:
float masking = 1. - smoothstep(sdfMax/2., sdfMax, abs(sdf));
if(false) {
for(int i=2; i>=0; i--) { bool rev = false;
//for(int i=0; i<3; i++) { bool rev = true;
p += rdSegmentRev(p, r, P[i], P[i+1], rev);
}
}
// JAGGED CASE (as rdrift superposition - symmetric but nonphysical)
if(true) {
bool rev = false;
vec2 u = vec2(0.);
u += rdSegmentRev(p, r, P[0], P[1], rev);
u += rdSegmentRev(p, r, P[1], P[2], rev);
u += rdSegmentRev(p, r, P[2], P[0], rev);
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).
vec2 rdSegmentKernel(vec2 p, float L) {
vec2 a = vec2(0., 0.);
vec2 b = vec2(L,0.);
float r = 1.;
// REVERSE COMB SETUP: b→a
vec2 rdir = a - b;
// ADVECT PARTICLE at p AS CYLINDER MOVES FROM b to a.
vec2 p0 = p; // Particle init (USED TO GET RevDrift: p-p0)
vec2 c = b; // Cylinder init
vec2 vc = normalize(rdir); // Cylinder unit velocity
// IF REMOVE TINE, INFLATE COORDS:
//p = applyTineInsertionMap(p, R ); // BODY FRAME
p = b + mxwApplyTineInsertionMap(p-b, 1.); // WORLD FRAME
// Reverse-advect p as cylinder moves from b→a:
vec2 q = p - c; // relative particle position
float dt = L;
vec2 vp = cylFlowVelWorldArbDir(q,vc); // particle velocity
// FWD EULER:
p += vp * dt; // position update
c += vc * dt; // advance cylinder
p = a + mxwApplyTineRemovalMap (p-a, 1.); // unit tine
vec2 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.
*/
vec2 rdTriWave(vec2 p, float r, vec2 dir, float L, float A, bool broken, bool dashed)
{
vec2 p0 = p; // initial position for reverse drift calc (p-p0)
dir = normalize(broken ? -dir : dir); // change direction so reversed rdSegments have same flow trend
// Project p onto line --> c and get coords for p as (x,y)
vec2 c = dot(p,dir)*dir;
float x = dot(c,dir);
float y = length(p - c); // unsigned dist to line
float segR = sqrt(A*A + L*L/16.);
float R = r * RDSEGMENT_MASK_R_FACTOR + segR;// Bounding radius of point w.r.t. qi
float H = L*0.5;
// 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.0 || H <= 0.0) return vec2(0.0);
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
vec2 vd = 0.5 *H *dir;
vec2 vp = A*rot90(dir);
for(int i=iR; i>=iL; i--) {// rev sweep
float si = (abs(i)%2==0) ? -1. : +1.; // -1even | +1odd
if(!dashed || (dashed && si<0.)) {
vec2 qi = dir*H*float(i);
vec2 ai = qi - vd + si*vp;
vec2 bi = qi + vd - si*vp;
p += rdSegmentRev(p, r, ai, bi, broken);
}
}
return p - p0; // reverse drift displacement
}
vec2 rdTriWaveComb(vec2 p, float r, vec2 dir, float L, float A, bool broken, bool dashed, float combGap)
{
dir = normalize(dir);
vec2 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.0) return vec2(0.0); // 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)
vec2 rdf = vec2(0.);
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
vec2 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 GLSL. Released under the MIT license.