← geometric-flows

Lattice flow (plane)

A planar lattice flowing under geodesic and horocyclic flows, tracing the orbit of its shortest vectors.

Source

#include "../shared/2d_utils.glsl"
#include "../shared/coord_utils.glsl"

bool showTrails = false;
float trailFalloff = 1.; // Length of fading trails, 1 gives no tails.
bool showDots = true;
bool showBasis = false;
bool secondLattice = false;
float traceLevel = 0.3;
float backFade = .999;
bool chBasis = true;

vec4 dotCol = vec4(.8,.05,.04,1.);
vec4 dotCol2 = vec4(.05,.04,.55,1.);
vec4 backCol = vec4(1.,1.,1.,0);
vec4 traceCol = vec4(0.,0.,0.,1.);
vec4 basisCol = vec4(.651,.669,.95,1.);

float sc = 14.;

// Buffer for trailing dots

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{    
    vec2 p = coords(fragCoord, iResolution, uScale, iResolution.xy * 0.5);
    
	vec2 uv = fragCoord / iResolution.xy;
    
    vec4 col = vec4(blur(dotTrails,uv,5.*iResolution.xy));
    if(showTrails) {
    	col = mix(col,vec4(backCol.xyz,0),iResolution.y/640.*trailFalloff/length(backCol-col));
    } else {
        col = vec4(backCol.xyz,0);
    }
    
    mat2 pLattice = mat2(uLattice);

    if(showDots) 
    { 
        col = drGrid(col, p, pLattice, 0.2, dotCol, 0.15); 
    }
    
    if(showBasis)
    {    
    	col = drPt(col, p, uLattice, .4, basisCol, .75);
    	col = drPt(col, p, vec4(uLattice.zw,0,0), .4, basisCol, .75);
    }

    fragColor = col;
}

// Buffer for Background and complete path record

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{   
    vec2 p = coords(fragCoord, iResolution, uScale, iResolution.xy * 0.5);
	vec2 uv = fragCoord / iResolution.xy;

    vec4 col = mix(backCol,vec4(texture(gridTrace,uv).xyz,1),backFade);

    col = drGrid(col, p, mat2(uLattice), 0.075, mix(col,traceCol,traceLevel), 2.);
    
    if(iFrame == 0 || uv.x > 1. || uv.y <0. || uv.y > 1. || uv.y < 0. || iMouse.z > 0.) 
    { col = backCol; }
    
    fragColor = col;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
	vec2 uv = fragCoord.xy / iResolution.xy;
    
    vec4 col = vec4(texture(gridTrace,uv).xyz,1);
    
    vec4 traced = vec4(texture(dotTrails,uv));
    
    col = vec4(mix(col, traced, traced.w).xyz,1.);

    // Output to screen
    fragColor = col;
}

//import * as math from 'mathjs';
import { matrix, dot, add, subtract, multiply, norm, hypot, subset, index, inv} from 'mathjs';
import * as cfUtils from '../shared/cf_utils.js';

let b1 = [1, 1.618033988749895];
let b2 = [-1.618033988749895, 1];

let p1 = null, p2 = null, hasReturnedFar = false;

let latticeSet = false;

function cDiv(a, b) {
  const d = b[0]*b[0] + b[1]*b[1];
  return [(a[0]*b[0] + a[1]*b[1]) / d, (a[1]*b[0] - a[0]*b[1]) / d];
}

function hpPt(b1, b2) {
  let pt = hypot(...b1) > hypot(...b2) ? cDiv(b1, b2) : cDiv(b2, b1);
  if (pt[1] < 0) pt = [-pt[0], -pt[1]];
  return pt;
}

function checkCF(flow) {
  console.log('checkPeriodSnap called, flow:', flow, 'p1:', p1);
  if (p1 === null) { return; }
  if (flow !== 0) { return; }

  

  const dist1 = norm(subtract(b1,p1));
  const dist2 = norm(subtract(b2,p2));

  const dist = Math.sqrt(dist1*dist1+dist2*dist2);
  

 if (dist < .01) {
    b1 = [...p1];
    b2 = [...p2];
    hasReturnedFar = false;
  }
}

// Input: basis as [[a, b], [c, d]] (rows are basis vectors)
// Output: reduced basis with b1 = shortest vector
function lagrangeGauss(b1, b2) {
  b1 = matrix(b1);
  b2 = matrix(b2);
  const normSq = v => dot(v, v);

  if (normSq(b1) > normSq(b2)) { [b1, b2] = [b2, b1]; }

  while (true) {
    const m = Math.round(dot(b1, b2) / normSq(b1));
    b2 = subtract(b2, multiply(m, b1));
    if (normSq(b2) >= normSq(b1)) break;
    [b1, b2] = [b2, b1];
  }

  if(subset(b1,index(1))<0) { b1 = multiply(-1,b1);}
  if(subset(b2,index(1))<0) { b2 = multiply(-1,b2);}

  return [b1.toArray(), b2.toArray()];
}

// Constrained Gauss–Lagrange reduction.
//   Keeps b1 in upper right quadrant
//   Keeps b2 in upper left quadrant
function coneReduce(b1, b2) {
  // Accept plain arrays or mathjs matrices (onFrame passes multiply results).
  b1 = (b1.toArray ? b1.toArray() : b1).slice();
  b2 = (b2.toArray ? b2.toArray() : b2).slice();
  const clamp = (k, lo, hi) => Math.max(lo, Math.min(hi, k));

  let changed = true;
  while (changed) {
    changed = false;

    // shorten b2 along b1, keeping b2 in QII (x <= 0, y >= 0)
    {
      const q = dot(b1, b1);
      if (q > 0) {
        let k = Math.round(-dot(b1, b2) / q);          // unconstrained minimiser
        const hi = b1[0] > 0 ? Math.floor(-b2[0] / b1[0]) :  Infinity; // x <= 0
        const lo = b1[1] > 0 ? Math.ceil (-b2[1] / b1[1]) : -Infinity; // y >= 0
        k = clamp(k, lo, hi);
        if (k !== 0) { b2 = add(b2, multiply(k, b1)); changed = true; }
      }
    }

    // shorten b1 along b2, keeping b1 in QI (x >= 0, y >= 0)
    {
      const q = dot(b2, b2);
      if (q > 0) {
        let k = Math.round(-dot(b1, b2) / q);
        const hi = b2[0] < 0 ? Math.floor(-b1[0] / b2[0]) :  Infinity; // x >= 0
        const lo = b2[1] > 0 ? Math.ceil (-b1[1] / b2[1]) : -Infinity; // y >= 0
        k = clamp(k, lo, hi);
        if (k !== 0) { b1 = add(b1, multiply(k, b2)); changed = true; }
      }
    }
  }
  return [b1, b2]; // b1/b2 are plain arrays (add/multiply preserve arrays)
}

const fps = 60;

// Geodesic (hyperbolic diagonal) flow matrix — stretches one axis, shrinks the other
function gFlow(dTime) {
  return matrix([
    [Math.exp(dTime / (3 * fps)), 0],
    [0, Math.exp(-dTime / (3 * fps))]
  ]);
}

// Horocyclic (shear) flow matrix (vertical)
function hFlow(dTime) {
  return matrix([
    [1, 0],
    [dTime / fps, 1]
  ]);
}

// Horocyclic (shear) flow matrix (horizontal)
function hFlowH(dTime) {
  return matrix([
    [1, dTime / fps],
    [0, 1]
  ]);
}

export function onSetLattice(inputs, handle) {
  const alpha = cfUtils.quadCFStrings(inputs.cfPrefix ?? '', inputs.cfPeriod ?? '0');
  const beta = cfUtils.quadCFStrings('', cfUtils.reverseCFString(inputs.cfPeriod) ?? '0');
  b1 = [1, beta];
  b2 = [-alpha, 1];

  [b1,b2] = coneReduce(b1,b2);
  
  

  //handle.setUniform('uLattice', [1, beta, -alpha, 1]);
  latticeSet = true;
}

export function setup(engine) {
  b1 = [1, 1.618033988749895];
  b2 = [-1.618033988749895, 1];
}

export function onFrame(engine, time, deltaTime, frame) {
  const flow = engine.getUniformValue('uFlow');
  const speed = engine.getUniformValue('uSpeed');

  let m;
  switch (flow) {
    case 1:  m = hFlow(speed); break;
    case 2:  m = hFlowH(speed); break;
    default: m = gFlow(speed * 2.);
  }

  const mb1 = multiply(m, b1);
  const mb2 = multiply(m, b2);

  [b1, b2] = coneReduce(mb1, mb2);

  checkCF(flow);

  // Live CF of the flow's forward endpoint — the future itinerary. This endpoint
  // ratio is invariant under the flow (the per-frame scale factors cancel), so it
  // stays stable and only advances one CF digit per coneReduce reduction (the
  // Gauss map). The invariant ratio depends on flow type/direction:
  //   geodesic (0): x-ratio forward when speed>=0, else y-ratio (reversed flow)
  //   horocyclic vertical (1):   x-ratio (parabolic fixed endpoint)
  //   horocyclic horizontal (2): y-ratio
  let endNum, endDen;
  if (flow === 2 || (flow === 0 && speed < 0)) {
    endNum = b2[1]; endDen = b1[1];   // y-ratio
  } else {
    endNum = b2[0]; endDen = b1[0];   // x-ratio
  }
  if (Math.abs(endDen) > 1e-9) {
    engine.emit('cf', cfUtils.cfExpansion(Math.abs(endNum / endDen), 12).join(', '));
  }

  // Inverse of this frame's smooth flow (pre-reduction), flattened column-major
  // so GLSL `mat2(uFlowInv) * p` walks a point back one flow step. Used by the
  // stretch-warp trails; shaders that don't read it simply ignore it.
  const mi = inv(m);
  engine.setUniformValue('uFlowInv', [mi.get([0, 0]), mi.get([1, 0]), mi.get([0, 1]), mi.get([1, 1])]);

  engine.setUniformValue('uLattice', [b1[0], b1[1], b2[0], b2[1]]);
  engine.setUniformValue('uNewLattice', latticeSet);
  latticeSet = false;
}