#!/usr/bin/env python3
"""A morel at night, with three visitors. Morel's first picture, Nov 9 2026.
Pure standard library: every pixel is computed, then written as a PNG by hand."""
import math, struct, zlib

W, H, SS = 600, 800, 2
CX = 300.0

def mix(a, b, t):
    t = max(0.0, min(1.0, t))
    return tuple(a[i] + (b[i] - a[i]) * t for i in range(3))

def smooth(e0, e1, x):
    t = max(0.0, min(1.0, (x - e0) / (e1 - e0)))
    return t * t * (3 - 2 * t)

def hash2(i, j):
    n = (i * 374761393 + j * 668265263) & 0xFFFFFFFF
    n = ((n ^ (n >> 13)) * 1274126177) & 0xFFFFFFFF
    return ((n ^ (n >> 16)) & 0xFFFF) / 65535.0

CAP_TOP, CAP_BOT = 120.0, 500.0
def cap_halfwidth(y):
    t = 1 - (y - CAP_TOP) / (CAP_BOT - CAP_TOP)   # pointed at the top, widest low down
    if t <= 0 or t >= 1:
        return 0.0
    return 128.0 * math.sin(math.pi * t ** 0.72) ** 0.72

def hexcell(px, py, r):
    """Nearest hexagon centre on a grid of radius r: returns (cell id, edge distance 0..1)."""
    sx, sy = r * 1.5, r * math.sqrt(3)
    best = None
    for off in (0, 1):
        ci = round((px - off * sx) / (2 * sx)) * 2 + off
        cj = round((py - (off * sy / 2)) / sy)
        qx = abs(px - ci * sx)
        qy = abs(py - (cj * sy + (off * sy / 2)))
        d = max(qx * 0.5 + qy * (math.sqrt(3) / 2), qx) / (r * math.sqrt(3) / 2 * 1.1547)
        if best is None or d < best[1]:
            best = ((ci, cj), d, (px - ci * sx, py - (cj * sy + off * sy / 2)))
    return best

BEES = [(455.0, 250.0, 0.35), (130.0, 395.0, -0.5), (480.0, 545.0, 0.9)]  # x, y, tilt

def bee(x, y, col):
    for bx, by, tilt in BEES:
        dx, dy = x - bx, y - by
        if abs(dx) > 30 or abs(dy) > 30:
            continue
        c, s = math.cos(tilt), math.sin(tilt)
        u, v = dx * c + dy * s, -dx * s + dy * c
        # wings: two pale ovals above the body
        for wx in (-5.0, 6.0):
            wd = ((u - wx) / 8.0) ** 2 + ((v + 11.0) / 11.0) ** 2
            if wd < 1:
                col = mix(col, (225, 235, 245), 0.55 * (1 - wd * 0.5))
        bd = (u / 14.0) ** 2 + (v / 9.0) ** 2
        if bd < 1:
            stripe = math.sin(u * 0.62 + 0.6) > 0.25
            body = (38, 28, 22) if stripe else (250, 196, 40)
            if u > 9:
                body = (38, 28, 22)   # head
            col = mix(body, (20, 14, 10), bd ** 3 * 0.6)
    return col

def shade(x, y):
    # night sky
    col = mix((9, 14, 34), (24, 44, 62), y / H)
    glow = math.exp(-(((x - CX) / 260) ** 2 + ((y - 330) / 300) ** 2))
    col = mix(col, (60, 80, 92), glow * 0.45)
    gi, gj = int(x // 23), int(y // 23)
    h = hash2(gi, gj)
    if h > 0.93 and y < 600:
        sx, sy = gi * 23 + 4 + hash2(gj, gi) * 15, gj * 23 + 4 + hash2(gi + 7, gj) * 15
        d = math.hypot(x - sx, y - sy)
        col = mix(col, (235, 235, 215), max(0, 1 - d / (0.8 + h * 1.2)) * 0.9)

    # ground
    gy = 668 + 12 * math.sin(x / 85.0) + 5 * math.sin(x / 23.0 + 1.3)
    if y > gy:
        col = mix((30, 52, 38), (12, 22, 20), (y - gy) / 130)
        if hash2(int(x // 3), int(y // 5)) > 0.8:
            col = mix(col, (52, 82, 50), 0.5)

    # stem
    if 440 < y < 705:
        t = (y - 440) / 265
        hw = 46 + 16 * t + 5 * math.sin(t * 6.0)
        if y > 680:
            hw *= math.sqrt(max(0.0, 1 - ((y - 680) / 25) ** 2))
        u = (x - CX - 4 * math.sin(t * 2.2)) / hw if hw > 0 else 9
        if abs(u) < 1:
            base = mix((238, 226, 196), (190, 172, 138), abs(u) ** 2.2)
            base = mix(base, (252, 244, 220), max(0, 0.5 - abs(u + 0.35)) * 0.7)
            if hash2(int(x // 2), int(y // 7)) > 0.86:
                base = mix(base, (205, 190, 158), 0.35)
            base = mix(base, (120, 104, 82), smooth(0.55, 0.0, (y - 440) / 70) * 0.55)  # shadow under the cap
            edge = smooth(1.0, 0.94, abs(u))
            col = mix(col, base, edge)

    # cap
    hw = cap_halfwidth(y)
    if hw > 0:
        u = (x - CX) / hw
        if abs(u) < 1:
            wrap = math.asin(max(-1, min(1, u))) * hw * 1.05      # bend the grid around the cap
            cell, d, (qx, qy) = hexcell(wrap + 3 * math.sin(y / 31.0), (y - CAP_TOP) * 0.66, 21.0)
            r = hash2(*cell)
            ridge = smooth(0.70 - r * 0.08, 0.92, d)
            pit = mix((58, 36, 22), (104, 68, 36), 0.25 + 0.75 * r)
            pit = mix(pit, (30, 18, 12), smooth(6, -14, qy) * 0.75)   # each pit is darker at its top
            rid = mix((214, 170, 96), (238, 204, 130), r)
            c = mix(pit, rid, ridge)
            c = mix(c, (255, 236, 170), ridge * max(0, 0.6 - abs(u + 0.4)) * 0.5)   # moonlit side
            c = mix(c, (40, 26, 18), smooth(0.55, 1.0, abs(u)) * 0.7)                # falls away at the edges
            col = mix(col, c, smooth(1.0, 0.95, abs(u)))
    return bee(x, y, col)

def render():
    rows = []
    inv = 1.0 / (SS * SS)
    for py in range(H):
        row = bytearray([0])
        for px in range(W):
            r = g = b = 0.0
            for sy in range(SS):
                for sx in range(SS):
                    c = shade(px + (sx + 0.5) / SS, py + (sy + 0.5) / SS)
                    r += c[0]; g += c[1]; b += c[2]
            row += bytes((int(r * inv), int(g * inv), int(b * inv)))
        rows.append(bytes(row))
    return b"".join(rows)

def png(path, raw):
    def chunk(tag, data):
        return struct.pack(">I", len(data)) + tag + data + struct.pack(">I", zlib.crc32(tag + data) & 0xFFFFFFFF)
    with open(path, "wb") as f:
        f.write(b"\x89PNG\r\n\x1a\n" + chunk(b"IHDR", struct.pack(">IIBBBBB", W, H, 8, 2, 0, 0, 0))
                + chunk(b"IDAT", zlib.compress(raw, 9)) + chunk(b"IEND", b""))

if __name__ == "__main__":
    png("morel-at-night.png", render())
    print("wrote morel-at-night.png")