// rd.c — a Gray-Scott reaction-diffusion sim that blooms in the terminal.
//
//   cc -O2 -o rd rd.c -lm && ./rd
//
// Two chemicals, u and v, sit on a toroidal grid. u is everywhere; a seed of v
// is dropped in the middle. Every step they diffuse and react:
//
//   u' = u + Du*Lap(u) - u*v*v + F*(1-u)
//   v' = v + Dv*Lap(v) + u*v*v - (F+k)*v
//
// That tiny feedback loop — v eats u to make more v, u trickles back in — is
// enough to grow the Turing patterns biology uses for spots, stripes and coral.
// Nothing here is drawn; the shapes are an emergent property of the numbers.
//
// Same trick as the path tracer next door: each terminal cell is two pixels
// stacked with the upper-half block '▀', foreground on top, background below.
// Same synthwave palette, too — violet field, magenta reaction, cyan on the
// fronts where the chemistry is moving fastest.

#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// ------------------------------------------------------------------------- rng
// xorshift64 — deterministic, so a given seed always blooms the same way.
static uint64_t rng_state = 0x2545f4914f6cdd1dULL;
static double rnd(void) {
    uint64_t x = rng_state;
    x ^= x << 13; x ^= x >> 7; x ^= x << 17;
    rng_state = x;
    return (x >> 11) * (1.0 / 9007199254740992.0); // [0,1)
}

// ------------------------------------------------------------------ vector type
typedef struct { double x, y, z; } V;
static V    cv(double x, double y, double z) { return (V){x, y, z}; }
static V    add(V a, V b)      { return cv(a.x+b.x, a.y+b.y, a.z+b.z); }
static V    mul(V a, double s) { return cv(a.x*s, a.y*s, a.z*s); }
static V    mix(V a, V b, double t) { return add(mul(a, 1-t), mul(b, t)); }

// ------------------------------------------------------------------ the chemistry
// Gray-Scott coefficients. Du/Dv are diffusion rates; F (feed) and k (kill) pick
// the regime. These defaults sit in the "coral" pocket: the seed branches out
// in growing fronts that fill the field with reef-like filaments.
static double Du = 1.0, Dv = 0.5;
static double F = 0.0545, k = 0.0620;

// 9-point Laplacian weights (orthogonal 0.2, diagonal 0.05, centre -1; sums to 0).
static void step(double *u, double *v, double *un, double *vn, int W, int H) {
    for (int y = 0; y < H; y++) {
        int yn = (y - 1 + H) % H, ys = (y + 1) % H;
        for (int x = 0; x < W; x++) {
            int xw = (x - 1 + W) % W, xe = (x + 1) % W;
            int i = y*W + x;
            // gather the 8 neighbours once, reuse for both fields
            int n = yn*W+x,  s = ys*W+x,  w = y*W+xw, e = y*W+xe;
            int nw = yn*W+xw, ne = yn*W+xe, sw = ys*W+xw, se = ys*W+xe;
            double lu = 0.2*(u[n]+u[s]+u[w]+u[e])
                      + 0.05*(u[nw]+u[ne]+u[sw]+u[se]) - u[i];
            double lv = 0.2*(v[n]+v[s]+v[w]+v[e])
                      + 0.05*(v[nw]+v[ne]+v[sw]+v[se]) - v[i];
            double uvv = u[i]*v[i]*v[i];
            un[i] = u[i] + Du*lu - uvv + F*(1.0 - u[i]);
            vn[i] = v[i] + Dv*lv + uvv - (F + k)*v[i];
        }
    }
}

// Seed: u=1 everywhere, then drop a few noisy squares of v near the centre.
static void seed(double *u, double *v, int W, int H) {
    for (int i = 0; i < W*H; i++) { u[i] = 1.0; v[i] = 0.0; }
    int cx = W/2, cy = H/2, r = (W < H ? W : H) / 14 + 2;
    for (int y = cy-r; y <= cy+r; y++)
        for (int x = cx-r; x <= cx+r; x++) {
            int i = ((y+H)%H)*W + ((x+W)%W);
            u[i] = 0.50 + 0.05*(rnd()-0.5);
            v[i] = 0.25 + 0.05*(rnd()-0.5);
        }
}

// ----------------------------------------------------------------------- palette
// Map a concentration (and the local front strength) to a synthwave colour:
// indigo where nothing's happening, violet -> hot magenta as v rises, and a
// cyan glow added on the reaction fronts where v changes sharply.
static V shade(double val, double front) {
    double t = val / 0.34; if (t < 0) t = 0; if (t > 1) t = 1;
    V bg  = cv(0.03, 0.01, 0.09);   // empty substrate
    V lo  = cv(0.30, 0.05, 0.55);   // violet
    V hi  = cv(1.05, 0.22, 0.62);   // hot magenta
    V col = t < 0.55 ? mix(bg, lo, t / 0.55)
                     : mix(lo, hi, (t - 0.55) / 0.45);
    double g = front * 7.0; if (g > 1) g = 1;
    col = add(col, mul(cv(0.10, 1.00, 1.30), g * 0.9)); // cyan front glow
    return col;
}

// linear -> sRGB-ish byte, matching the path tracer's tonemap so the two
// exhibits sit at the same exposure.
static int tone(double c) {
    if (c < 0) c = 0;
    c = pow(c / (c + 1.0), 1.0/2.2);
    int x = (int)(c*255 + 0.5);
    return x < 0 ? 0 : x > 255 ? 255 : x;
}

// Build an RGB framebuffer from the v field. front = |v - mean(4-neighbours)|.
static void colorize(double *v, V *fb, int W, int H) {
    for (int y = 0; y < H; y++)
        for (int x = 0; x < W; x++) {
            int i = y*W + x;
            int n = ((y-1+H)%H)*W+x, s = ((y+1)%H)*W+x;
            int w = y*W+((x-1+W)%W),  e = y*W+((x+1)%W);
            double front = fabs(v[i] - 0.25*(v[n]+v[s]+v[w]+v[e]));
            fb[i] = shade(v[i], front);
        }
}

static void write_ppm(const char *path, V *fb, int W, int H) {
    FILE *fp = fopen(path, "wb");
    fprintf(fp, "P6\n%d %d\n255\n", W, H);
    for (int i = 0; i < W*H; i++) {
        unsigned char px[3] = { tone(fb[i].x), tone(fb[i].y), tone(fb[i].z) };
        fwrite(px, 1, 3, fp);
    }
    fclose(fp);
}

int main(int argc, char **argv) {
    int W = 80, H = 48;          // grid / framebuffer size
    int SPF = 14;                // simulation steps per displayed frame
    int FRAMES = 320;
    const char *snapshot = NULL; // --shot: settle, then write one PPM
    const char *gifpre = NULL;   // --gif: dump every frame as a PPM
    if (argc > 1) SPF = atoi(argv[1]);
    if (argc > 2) FRAMES = atoi(argv[2]);
    if (argc > 3 && !strcmp(argv[3], "--shot")) {
        W = atoi(argv[4]); H = atoi(argv[5]); snapshot = argv[6];
    }
    if (argc > 3 && !strcmp(argv[3], "--gif")) {
        W = atoi(argv[4]); H = atoi(argv[5]); gifpre = argv[6];
    }

    double *u  = malloc(sizeof(double)*W*H), *v  = malloc(sizeof(double)*W*H);
    double *un = malloc(sizeof(double)*W*H), *vn = malloc(sizeof(double)*W*H);
    V *fb = malloc(sizeof(V)*W*H);
    char *out = malloc(W*H*48 + 4096);
    seed(u, v, W, H);

    // A hidden still: run the chemistry to a rich pattern, then dump one frame.
    if (snapshot) {
        for (int s = 0; s < FRAMES*SPF; s++) {
            step(u, v, un, vn, W, H);
            double *tu=u; u=un; un=tu; double *tv=v; v=vn; vn=tv;
        }
        colorize(v, fb, W, H);
        write_ppm(snapshot, fb, W, H);
        return 0;
    }

    if (!gifpre) { printf("\033[2J\033[?25l"); fflush(stdout); }

    for (int f = 0; f < FRAMES; f++) {
        for (int s = 0; s < SPF; s++) {
            step(u, v, un, vn, W, H);
            double *tu=u; u=un; un=tu; double *tv=v; v=vn; vn=tv;
        }
        colorize(v, fb, W, H);

        if (gifpre) {
            char path[512];
            snprintf(path, sizeof path, "%s%03d.ppm", gifpre, f);
            write_ppm(path, fb, W, H);
            fprintf(stderr, "\rframe %d/%d", f+1, FRAMES);
            continue;
        }

        char *o = out;
        o += sprintf(o, "\033[H");
        for (int y = 0; y < H; y += 2) {
            for (int x = 0; x < W; x++) {
                V top = fb[y*W + x], bot = fb[(y+1)*W + x];
                o += sprintf(o, "\033[38;2;%d;%d;%dm\033[48;2;%d;%d;%dm▀",
                             tone(top.x), tone(top.y), tone(top.z),
                             tone(bot.x), tone(bot.y), tone(bot.z));
            }
            o += sprintf(o, "\033[0m\n");
        }
        fwrite(out, 1, o - out, stdout);
        fflush(stdout);
    }

    if (gifpre) fprintf(stderr, "\n");
    else printf("\033[0m\033[?25h\n");
    free(u); free(v); free(un); free(vn); free(fb); free(out);
    return 0;
}