Hey, it kinda works! Single-threaded

diff --git a/main.c b/main.c
new file mode 100755 (executable)
index 0000000..58c65fb
--- /dev/null
+++ b/main.c
@@ -0,0 +1,191 @@
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#define GRID_X 256
+#define GRID_Y 256
+
+typedef struct
+{
+    float v[3][3];
+
+} Mat3;
+
+typedef struct
+{
+    float a;
+    float b;
+} FVec2;
+
+typedef struct
+{
+    float diff_a;
+    float diff_b;
+    float feed;
+    float kill;
+    float delta_t;
+} RD_Opts;
+
+FVec2 grid[GRID_X][GRID_Y];
+FVec2 grid_prime[GRID_X][GRID_Y];
+
+const Mat3 laplacian_kernel =
+{
+    // oh shit
+    // Sum of all members must balance to 0.0 (or outer members have to do 1.0)
+    // it looks like any modifications should be on the order of 0.01 or less
+    .v = {
+        {0.05f, 0.2f, 0.05f},
+        {0.20, -1.0f, 0.20},
+        {0.05f, 0.2f, 0.05f}
+    }
+};
+
+float convolve(Mat3 kernel, Mat3 source)
+{
+    float result = 0.0f;
+
+    // Convolve kernel at grid[x,y]
+    for (int ix = 0; ix < 3; ix++)
+    {
+        for (int iy = 0; iy < 3; iy++)
+        {
+            result += (kernel.v[ix][iy] * source.v[ix][iy]);
+        }
+    }
+
+    return result;
+}
+
+float rd_a_prime(RD_Opts opts, int x, int y, Mat3 kernel, float A, float B)
+{
+    float a_prime = 1.0f;
+    // Use species A in the convolution, b_prime will need B
+    // For now we won't iterate over edge rows and columns to avoid special case vomit
+    Mat3 source =
+    {
+        .v =
+        {
+            {grid[x-1][y-1].a,  grid[x][y-1].a, grid[x+1][y-1].a},
+            {grid[x-1][y+0].a,  grid[x][y+0].a, grid[x+1][y+0].a},
+            {grid[x-1][y+1].a,  grid[x][y+1].a, grid[x+1][y+1].a},
+        }
+    };
+
+    a_prime = A + (opts.diff_a*(convolve(kernel, source)) - (A*B*B) + opts.feed*(1.0f-A)) * opts.delta_t;
+    if(isnan(a_prime))
+    {
+        int bp = 0;
+        printf("nan!\n");
+    }
+    return a_prime;
+}
+
+float rd_b_prime(RD_Opts opts, int x, int y, Mat3 kernel, float A, float B)
+{
+    float b_prime = 1.0f;
+    // Use species A in the convolution, b_prime will need B
+    Mat3 source =
+    {
+        .v =
+        {
+            {grid[x-1][y-1].b,  grid[x][y-1].b, grid[x+1][y-1].b},
+            {grid[x-1][y+0].b,  grid[x][y+0].b, grid[x+1][y+0].b},
+            {grid[x-1][y+1].b,  grid[x][y+1].b, grid[x+1][y+1].b},
+        }
+    };
+
+    b_prime = B + (opts.diff_b*(convolve(kernel, source)) + A*(B*B) - (opts.kill + opts.feed)*B) * opts.delta_t;
+
+    if(isnan(b_prime))
+    {
+        int bp = 0;
+        printf("nan!\n");
+    }
+    return b_prime;
+}
+
+int main(int argc, char** argv)
+{
+    char chars[6] = {' ', ' ', ' ', '+', '#', '@'};
+    const int max_chars = sizeof(chars) / sizeof(chars[0]) - 1;
+    // Connectivity increases as time passes (iterations or time delta)?
+    // Seems we like <1.0 time deltas and higher iterations
+
+    RD_Opts puffer_opts =
+    {
+        .diff_a = 1.0f,
+        .diff_b = 0.5f,
+        .feed = 0.055f,
+        .kill = 0.062f,
+        .delta_t = 0.6f,
+    };
+
+    RD_Opts coral_opts =
+    {
+        .diff_a = 0.7f,
+        .diff_b = 0.3f,
+        .feed = 0.0545f,
+        .kill = 0.062f,
+        .delta_t = 1.0f,
+    };
+
+               RD_Opts opts = coral_opts;
+
+
+    for (int x = 0; x < GRID_X; x++)
+    {
+        for (int y = 0; y < GRID_Y; y++)
+        {
+            grid[x][y].a = 1.0f;
+            grid[x][y].b = 0.0f;
+        }
+    }
+
+    grid[GRID_X/2][GRID_Y/2].b = 1.0f;
+    grid[GRID_X/4][GRID_Y/2].b = 1.0f;
+    grid[(GRID_X/6) * 5][4 * (GRID_Y/9)].b = 1.0f;
+    grid[GRID_X/4 * 3][GRID_Y/3].b = 1.0f;
+    grid[GRID_X-1][GRID_Y/5].b = 1.0f;
+
+    const int max_iterations = 1 * 1e4;
+    for (int iterations = 0; iterations < max_iterations; iterations++)
+    {
+                       for (int x = 1; x < GRID_X-1; x++)
+                       {
+                               for (int y = 1; y < GRID_Y-1; y++)
+                               {
+                                       FVec2 each = grid[x][y];
+                                       grid_prime[x][y].a = rd_a_prime(opts, x, y, laplacian_kernel, each.a, each.b);
+                                       grid_prime[x][y].b = rd_b_prime(opts, x, y, laplacian_kernel, each.a, each.b);
+                               }
+                       }
+                       if (iterations % (int)1e3 == 0)
+                       {
+                               for (int x = 1; x < GRID_X-1; x++)
+                               {
+                                       for (int y = 1; y < GRID_Y-1; y++)
+                                       {
+                                               int idx = floor(grid[x][y].a * (max_chars));
+                                               printf("%c", chars[idx]);
+                                       }
+                                       printf("\n");
+                               }
+                       }
+                       memcpy(grid, grid_prime, sizeof(FVec2) * GRID_X * GRID_Y);
+    }
+
+    printf("Opts: {\nIterations: %d\nTime delta: %f\nDiffA: %f\nDiffB: %f\nFeed: %f\nKill: %f\n\n}\n", max_iterations, opts.delta_t, opts.diff_a, opts.diff_b, opts.feed, opts.kill);
+
+    for (int x = 1; x < GRID_X-1; x++)
+    {
+        for (int y = 1; y < GRID_Y-1; y++)
+        {
+            int idx = floor(grid[x][y].a * (max_chars));
+            printf("%c", chars[idx]);
+        }
+        printf("\n");
+    }
+
+    return 0;
+}
+