Canvas.setpenopacity(1);

const cells = 500; // min=250 max=750 step=50 Surface resolution
const pores = 3.8; // min=2.5 max=5 step=0.1 Pore frequency
const wall = 0.31; // min=0.2 max=0.6 step=0.01 Shell thickness
const twist = 0.2; // min=-0.5 max=0.5 step=0.05 Twist per unit height
const yaw = 30; // min=0 max=180 step=5 Horizontal view angle
const pitch = 38; // min=-60 max=60 step=5 Vertical view angle
const spacing = 0.045; // min=0.035 max=0.08 step=0.005 Engraving spacing

const turtle = new Turtle();
const {sin, cos, abs, max, min, floor, ceil, hypot} = Math;
const cosYaw = cos(yaw * Math.PI / 180);
const sinYaw = sin(yaw * Math.PI / 180);
const cosPitch = cos(pitch * Math.PI / 180);
const sinPitch = sin(pitch * Math.PI / 180);

const rowWidth = cells + 1;
const gridStep = 5.6 / cells;
const drawingScale = 32;
const samples = [];
const vertices = new Map();
const segments = [];
let contourPaths;
let nextPath = 0;

function toObjectSpace(x, y, z) {
    const u = x * cosYaw + z * sinYaw;
    const w = -x * sinYaw + z * cosYaw;
    return [u, y * cosPitch - w * sinPitch, y * sinPitch + w * cosPitch];
}

function sceneDistance(x, y, z) {
    const p = toObjectSpace(x, y, z);
    const angle = twist * p[1];
    const u = pores * (p[0] * cos(angle) - p[2] * sin(angle));
    const v = pores * p[1];
    const w = pores * (p[0] * sin(angle) + p[2] * cos(angle));

    // A gyroid repeats forever. The ellipsoid cuts out one piece of it.
    // Taking abs(g) keeps a thin wall on both sides of its zero surface.
    const g = sin(u) * cos(v) + sin(v) * cos(w) + sin(w) * cos(u);
    const shell = (abs(g) - wall) / (pores * 2.5);
    const ellipsoid = (hypot(p[0] / .92, p[1] / 1.13, p[2] / .92) - 2.35) * .8;
    return max(shell, ellipsoid);
}

function sampleSurface(x, y) {
    let z = 4.3;
    let hit = false;

    // This field only approximates distance, so take smaller steps to
    // reduce the chance of stepping straight through a thin wall.
    for (let i = 0; i < 260 && z > -4.3; i++) {
        const d = sceneDistance(x, y, z);
        if (d < .0008) {
            hit = true;
            break;
        }
        z -= max(.0005, d * .8);
    }
    if (!hit) return {x, y, hit: 0};

    const epsilon = .002;
    let nx = sceneDistance(x + epsilon, y, z) - sceneDistance(x - epsilon, y, z);
    let ny = sceneDistance(x, y + epsilon, z) - sceneDistance(x, y - epsilon, z);
    let nz = sceneDistance(x, y, z + epsilon) - sceneDistance(x, y, z - epsilon);
    const normalLength = hypot(nx, ny, nz);
    nx /= normalLength;
    ny /= normalLength;
    nz /= normalLength;

    // Probe just outside the surface. Nearby walls darken the recesses.
    let openness = 0;
    for (const [h, weight] of [[.03, .6], [.08, .25], [.18, .15]]) {
        const clearance = sceneDistance(x + nx * h, y + ny * h, z + nz * h) / h;
        openness += max(0, min(1, clearance)) * weight;
    }
    const light = max(.04, (-.5 * nx + .65 * ny + .65 * nz) / 1.047) * (.25 + .75 * openness);
    const p = toObjectSpace(x, y, z);

    // These planes live in object space, so their contours wrap around
    // the shell instead of looking like stripes laid over the picture.
    return {
        x, y, z, hit: 1, light,
        bands: .75 * p[1] + .25 * p[2],
        hatch: .6 * p[0] - .3 * p[1] + .5 * p[2]
    };
}

function vertexAt(point, family) {
    // Neighbouring triangles calculate the same endpoint independently.
    // Round the lookup key to join them despite tiny floating-point differences.
    const key = family + ':' + Math.round(point[0] * 1e6) + ',' + Math.round(point[1] * 1e6);
    if (!vertices.has(key)) vertices.set(key, {point, edges: []});
    return vertices.get(key);
}

function addSegment(a, b, family) {
    const start = vertexAt(a, family);
    const end = vertexAt(b, family);
    if (start === end) return;

    const edge = {start, end, used: false};
    start.edges.push(edge);
    end.edges.push(edge);
    segments.push(edge);
}

function traceContour(triangle, field, level, family, cutoff = 2) {
    const crossings = [];
    for (let i = 0; i < 3; i++) {
        const a = triangle[i];
        const b = triangle[(i + 1) % 3];
        if ((a[field] < level) === (b[field] < level)) continue;

        const t = (level - a[field]) / (b[field] - a[field]);
        crossings.push([
            a.x + (b.x - a.x) * t,
            a.y + (b.y - a.y) * t,
            field === 'hit' ? 0 : a.light + (b.light - a.light) * t
        ]);
    }
    if (crossings.length !== 2) return;

    // Extra hatch lines stop at the light threshold. Clipping here avoids
    // a whole triangle disappearing when just one corner is too bright.
    let [a, b] = crossings;
    if (a[2] >= cutoff && b[2] >= cutoff) return;
    if ((a[2] < cutoff) !== (b[2] < cutoff)) {
        const t = (cutoff - a[2]) / (b[2] - a[2]);
        const clipped = [a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1]), cutoff];
        if (a[2] >= cutoff) a = clipped;
        else b = clipped;
    }
    addSegment(a.slice(0, 2), b.slice(0, 2), family);
}

function traceBands(triangle, field, interval, offset, family, cutoff) {
    const low = min(...triangle.map(p => p[field]));
    const high = max(...triangle.map(p => p[field]));
    for (let k = ceil((low - offset) / interval); k <= floor((high - offset) / interval); k++) {
        traceContour(triangle, field, k * interval + offset, family, cutoff);
    }
}

function engraveTriangle(triangle) {
    traceContour(triangle, 'hit', .5, 'outline');
    if (triangle.some(p => !p.hit)) return;

    // A sudden depth change belongs to a different surface behind this one.
    // Don't bridge that gap with an engraving line.
    const depthJump = max(...triangle.map(p => p.z)) - min(...triangle.map(p => p.z));
    if (depthJump > .12) return;

    traceBands(triangle, 'bands', spacing, 0, 'main', 2);
    traceBands(triangle, 'bands', spacing, spacing / 2, 'shadow', .48);
    traceBands(triangle, 'hatch', spacing * .85, 0, 'cross', .24);
}

function joinContours() {
    const result = [];

    function follow(start, edge) {
        const path = [start.point];
        let at = start;
        while (edge && !edge.used) {
            edge.used = true;
            at = edge.start === at ? edge.end : edge.start;
            path.push(at.point);
            edge = at.edges.find(e => !e.used);
        }
        if (path.length > 3) result.push(path);
    }

    // Start with the open contours. Anything left afterwards is a loop.
    for (const vertex of vertices.values()) {
        if (vertex.edges.length === 2) continue;
        for (const edge of vertex.edges) {
            if (!edge.used) follow(vertex, edge);
        }
    }
    for (const edge of segments) {
        if (!edge.used) follow(edge.start, edge);
    }
    return result;
}

function drawContour(path) {
    // Two corner-cutting passes soften the sampling grid. Keep the ends
    // fixed so smoothing doesn't pull the lines away from the silhouette.
    let points = path;
    for (let pass = 0; pass < 2; pass++) {
        const smooth = [points[0]];
        for (let j = 1; j < points.length; j++) {
            const a = points[j - 1];
            const b = points[j];
            smooth.push([.75 * a[0] + .25 * b[0], .75 * a[1] + .25 * b[1]]);
            smooth.push([.25 * a[0] + .75 * b[0], .25 * a[1] + .75 * b[1]]);
        }
        smooth.push(points[points.length - 1]);
        points = smooth;
    }

    turtle.penup();
    turtle.goto(points[0][0] * drawingScale, -points[0][1] * drawingScale);
    turtle.pendown();
    for (let j = 1; j < points.length; j++) {
        turtle.goto(points[j][0] * drawingScale, -points[j][1] * drawingScale);
    }
}

function walk(i) {
    // Sample first, trace next, then draw. Doing one row at a time lets
    // Turtletoy yield between chunks of the expensive surface calculation.
    if (i < rowWidth) {
        for (let x = 0; x < rowWidth; x++) {
            samples.push(sampleSurface(-2.8 + x * gridStep, -2.8 + i * gridStep));
        }
    } else if (i < rowWidth + cells) {
        const y = i - rowWidth;
        for (let x = 0; x < cells; x++) {
            const a = samples[y * rowWidth + x];
            const b = samples[y * rowWidth + x + 1];
            const c = samples[(y + 1) * rowWidth + x + 1];
            const d = samples[(y + 1) * rowWidth + x];
            engraveTriangle([a, b, c]);
            engraveTriangle([a, c, d]);
        }
    } else {
        if (!contourPaths) contourPaths = joinContours();
        if (nextPath === contourPaths.length) return false;
        drawContour(contourPaths[nextPath++]);
    }
    return true;
}