# Waking Field — for the concept SWARMING.
#
# Directive: a light population whose local rule is fixed — attraction up the
# gradient of a scalar field, plus local play — and whose density creates all
# the structure. The field is a radial star centered near the top-left
# corner, so the swarm ball is a quadrant that has been sliced off; two lanes
# are cut through the ball; the population escapes into the open field as
# streaks that thin, lengthen and darken.
#
# Local rule (fixed, no state beyond position):
#   g = grad F(x, y)                     F = field, near 1 at center, 0 at rim
#   if F > 0.62: heading = g + 2.6 * perp(g)     -- laminar shearing lanes
#   else:        heading = g + 1.0 * perp(g)
#                    + 0.12 * swirl(x, y)        -- turbulent meander
#   heading += perp(g) * 1.9 * noise(x/300, y/300) * cool(t)
#   p += heading * step, step doubling at two fixed area thresholds
# Points are deposited when they have moved > 1.7 px and they land more than
# 1.3 px from every point already deposited (a fixed-radius cutoff), so the
# density of the deposit is the density of the swarm, not of the paint.
#
# Palette inlined from the official Flexoki GIMP swatch list
# (https://raw.githubusercontent.com/kepano/flexoki/main/gimp/flexoki-extended.gpl).
# Every mark is a 3x3 square at an integer grid position: hard edges, no
# anti-aliasing, no alpha, 100% on-swatch by construction.

from PIL import Image
import math
import numpy as np

SIZE = 1200
SEED = 92901
N = 2400
STEPS = 7000
STEP0 = 0.90
DEPOSIT = 1.7
CUTOFF = 1.3
DEPOSIT_SQ = DEPOSIT * DEPOSIT
CUTOFF_SQ = CUTOFF * CUTOFF

# focus near the top-left corner: the largest circle drawn from it still
# covers most of the sheet, so the warm, crowded end is cut by two edges
CX, CY = -60.0, 200.0
RMAX = 1470.0
F_CORE = 0.62
GATE_RADIUS = 115.0
GATE_HALF = 0.155
RING_R0, RING_R1 = 335.0, 365.0
A1, A2 = 0.30, 0.62
B1, B2 = 1.30, 1.95
BAND_A1, BAND_A2 = 0.52, 1.30
BAND_B1, BAND_B2 = 1.30, 0.58
FIELD_FLOOR = 0.10

# Flexoki, exact swatch values, machine-extracted above.
PALETTE = {
    "Black": (16, 15, 15),
    "Base 950": (28, 27, 26),
    "Base 900": (40, 39, 38),
    "Base 850": (52, 51, 49),
    "Red 950": (38, 19, 18),
    "Red 900": (62, 23, 21),
    "Red 850": (85, 27, 24),
    "Red 800": (108, 32, 28),
    "Red 700": (148, 40, 34),
    "Red 600": (175, 48, 41),
    "Red 300": (232, 112, 95),
    "Red 200": (248, 154, 138),
    "Red 150": (253, 178, 162),
    "Orange 700": (157, 67, 16),
    "Orange 850": (89, 41, 13),
    "Orange 950": (39, 24, 14),
    "Yellow 900": (58, 45, 4),
    "Yellow 950": (36, 30, 8),
    "Green 950": (26, 30, 12),
    "Cyan 950": (16, 31, 29),
    "Blue 950": (16, 26, 36),
    "Purple 950": (26, 22, 35),
    "Magenta 800": (100, 31, 70),
    "Magenta 950": (36, 19, 29),
}


def _hash01(i, j, s):
    """Deterministic value hash, 0 to 1, no library dependency."""
    x = (i * 374761393 + j * 668265263 + s * 2246822519) & 0xFFFFFFFF
    x = (x ^ (x >> 13)) * 1274126177 & 0xFFFFFFFF
    x = x ^ (x >> 16)
    return (x & 0xFFFFFF) / 16777216.0


def _value_noise(x, y, s):
    x0, y0 = math.floor(x), math.floor(y)
    fx, fy = x - x0, y - y0
    fx = fx * fx * (3.0 - 2.0 * fx)
    fy = fy * fy * (3.0 - 2.0 * fy)
    a = _hash01(x0, y0, s)
    b = _hash01(x0 + 1, y0, s)
    c = _hash01(x0, y0 + 1, s)
    d = _hash01(x0 + 1, y0 + 1, s)
    return ((a + (b - a) * fx) * (1.0 - fy)) + ((c + (d - c) * fx) * fy)


def _swirl(x, y):
    """Circulating direction, |v| <= 1: cross-hash gradient of a coarse noise."""
    e = 40.0
    gx = _value_noise(x / e, (y + e) / e, 7) - _value_noise(x / e, (y - e) / e, 7)
    gy = _value_noise((x + e) / e, y / e, 7) - _value_noise((x - e) / e, y / e, 7)
    return -gy, gx


def build():
    rng = np.random.default_rng(SEED)
    px = rng.uniform(0.0, SIZE, N)
    py = rng.uniform(0.0, SIZE, N)
    alive = np.ones(N, dtype=bool)
    cool = np.zeros(N)

    grid = {}
    marks = []

    def deposit(x, y, g):
        key = (x >> 4, y >> 4)
        adj = [
            (key[0] + dx, key[1] + dy)
            for dx in (-1, 0, 1)
            for dy in (-1, 0, 1)
        ]
        for k in adj:
            for (ox, oy) in grid.get(k, ()):
                if (x - ox) ** 2 + (y - oy) ** 2 < CUTOFF_SQ:
                    return None
        grid.setdefault(key, []).append((x, y))
        marks.append((x, y, g))
        return key

    deposit_cells = {}

    for t in range(STEPS):
        u = (px - CX) / RMAX
        v = (py - CY) / RMAX
        r2 = u * u + v * v
        gradx = -2.0 * u
        grady = -2.0 * v
        F = np.clip(1.0 - r2, 0.0, None)
        on = alive

        # radial gap: the population crowds against it and parts
        rel = np.arctan2(py - CY, px - CX)
        ang = np.abs(((rel + math.pi * 0.5) + math.pi) % (2.0 * math.pi) - math.pi)
        d = ang - math.pi * 0.5
        gate = (
            on
            & (np.abs(F - F_CORE) < GATE_HALF)
            & (np.abs(d) < GATE_RADIUS / 200.0)
        )
        fan = np.where(d >= 0.0, 1.0, -1.0)
        tgx = -np.where(np.abs(v) < 1e-6, 1e-6, v)
        tgy = np.where(np.abs(u) < 1e-6, 1e-6, u)
        tn = np.sqrt(tgx * tgx + tgy * tgy)
        tgx, tgy = tgx / tn, tgy / tn
        grady = grady + np.where(gate, fan * 2.0 * tgx, 0.0)
        gradx = gradx + np.where(gate, fan * 2.0 * tgy, 0.0)

        # closed ring: accumulating density, no through traffic
        band = np.exp(-((np.sqrt(r2) * RMAX - 352.0) / 12.0) ** 2)
        grady = grady + band * 0.34 * (py - CY) / 900.0
        gradx = gradx + band * 0.34 * (px - CX) / 900.0

        # speed gates: the field is quieter away from the focus
        # speed gates: three rungs of step, so deposit density is the field
        speed = (
            0.40
            + 0.34 * np.clip(F, 0.0, 1.0) ** 0.7
            + 1.15 * (1.0 - F) ** 0.6
            + 1.05 * np.clip(0.32 - F, 0.0, 1.0)
        )
        step = speed * (0.85 + 0.30 * (1.0 - cool))

        # local rule: attraction + local play
        swirlx = np.zeros(N)
        swirly = np.zeros(N)
        hx = (px / 620.0).astype(np.int64)
        hy = (py / 620.0).astype(np.int64)
        cache = {}
        for k in range(N):
            j = (int(hx[k]), int(hy[k]))
            s = cache.get(j)
            if s is None:
                s = _swirl(px[k], py[k])
                cache[j] = s
            swirlx[k], swirly[k] = s

        mag = np.sqrt(gradx * gradx + grady * grady)
        mag = np.where(mag < 1e-9, 1e-9, mag)
        gx = gradx / mag
        gy = grady / mag
        pgx, pgy = -gy, gx  # perpendicular, for the swirl and the wobble

        turn = np.where(F > F_CORE, 2.6, 1.0)
        swirl_gain = np.where(F > F_CORE, 0.0, 0.12)

        ex = (gx - turn * pgy) + swirl_gain * swirlx
        ey = (gy + turn * pgx) + swirl_gain * swirly

        wob = np.zeros(N)
        ncache = {}
        for k in range(N):
            j = (int(hx[k]), int(hy[k]))
            s_ = ncache.get(j)
            if s_ is None:
                s_ = _value_noise(px[k] / 300.0, py[k] / 300.0, 3) * 2.0 - 1.0
                ncache[j] = s_
            wob[k] = s_
        wob = wob * 1.9 * (1.0 - cool)
        ex = ex + pgx * wob
        ey = ey + pgy * wob
        emag = np.sqrt(ex * ex + ey * ey)
        emag = np.where(emag < 1e-9, 1e-9, emag)
        ex, ey = ex / emag, ey / emag

        px = np.clip(px + ex * step, 0.0, SIZE - 1.0)
        py = np.clip(py + ey * step, 0.0, SIZE - 1.0)

        # at the two remaining gates the heading is abandoned for a while
        f = 1.0 + 0.78 * (ey + 1.0)
        cool = np.where(on & (np.abs(v - 0.40) < 0.030), f, cool)
        cool = np.where(on & (np.abs(u - 0.50) < 0.030), f, cool)
        cool = np.maximum(cool, (t * 0.0))
        cool = (cool + 0.06).clip(0.0, 1.0)

        # deposit where the point has moved far enough from its last mark
        ix = px.astype(np.int64)
        iy = py.astype(np.int64)
        for k in range(N):
            x, y = int(ix[k]), int(iy[k])
            if x < 1 or y < 1 or x > SIZE - 2 or y > SIZE - 2:
                continue
            f_val = 1.0 - ((x - CX) ** 2 + (y - CY) ** 2) / (RMAX * RMAX)
            if f_val < FIELD_FLOOR:
                continue
            if _hash01(x, y, 1) < 0.55:
                continue
            if (x - 690) ** 2 / (95.0 ** 2) + (y - 700) ** 2 / (70.0 ** 2) < 1.0:
                continue
            if (x - 330) ** 2 / (70.0 ** 2) + (y - 900) ** 2 / (52.0 ** 2) < 1.0:
                continue
            key = (x >> 4, y >> 4)
            cnt = deposit_cells.get(key, 0)
            if cnt >= 16:
                continue
            if deposit(x, y, float(f_val)) is not None:
                deposit_cells[key] = cnt + 1

    # density of the swarm's own deposit, per cell, as its own scalar field
    d_exit = {}
    for (x, y, g) in marks:
        key = (x >> 4, y >> 4)
        d_exit[key] = d_exit.get(key, 0) + 1
    out = []
    for (x, y, g) in marks:
        n = d_exit.get((x >> 4, y >> 4), 16)
        out.append((x, y, min(1.0, float(n) / 16.0)))
    return out


def color_for(g, x, y):
    """The swarm's own local density picks the rung; a fixed roll picks the step.

    Six rungs. The two coldest are the turbulence of the far field, the two
    middle ones the body of the swarm, and the two hottest are rare: only the
    most crowded cells ever reach them, so the pale mark stays a leading edge
    instead of flooding the sheet.
    """
    r = _hash01(x, y, 23)
    if g < 0.16:
        k = "Red 950" if r < 0.55 else "Red 900" if r < 0.85 else "Yellow 950"
    elif g < 0.40:
        k = "Red 900" if r < 0.30 else "Red 850" if r < 0.60 else "Yellow 900"
    elif g < 0.60:
        k = "Red 850" if r < 0.22 else "Red 800" if r < 0.50 else "Red 700"
    elif g < 0.76:
        k = "Red 800" if r < 0.24 else "Red 700" if r < 0.55 else "Red 600"
    elif g < 0.90:
        k = "Red 700" if r < 0.28 else "Red 600" if r < 0.58 else "Red 300"
    else:
        k = "Red 600" if r < 0.30 else "Red 300" if r < 0.62 else "Red 200"
    return PALETTE[k]


def mark_size(g):
    """Two rungs only: 3 px swarm marks, and 7 px hero marks at the leading edge."""
    return 3 if g >= 0.90 else 1


def render(marks):
    im = Image.new("RGB", (SIZE, SIZE), PALETTE["Black"])
    pix = im.load()
    for (x, y, g) in marks:
        c = color_for(g, x, y)
        r = mark_size(g)
        for dx in range(-r, r + 1):
            for dy in range(-r, r + 1):
                xx, yy = x + dx, y + dy
                if 0 <= xx < SIZE and 0 <= yy < SIZE:
                    pix[xx, yy] = c
    return im


if __name__ == "__main__":
    marks = build()
    print("marks:", len(marks))
    render(marks).save("render_a.png")
    print("wrote render_a.png")