--- /dev/null
+#include "structs.h"
+
+FVec2 grid[GRID_X][GRID_Y];
+FVec2 grid_prime[GRID_X][GRID_Y];
+
+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;
+
+ 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;
+
+ return b_prime;
+}
+
--- /dev/null
+
+float convolve(Mat3 kernel, Mat3 source);
+float rd_a_prime(RD_Opts opts, int x, int y, Mat3 kernel, float A, float B);
+float rd_b_prime(RD_Opts opts, int x, int y, Mat3 kernel, float A, float B);
+
+extern FVec2 grid[GRID_X][GRID_Y];
+extern FVec2 grid_prime[GRID_X][GRID_Y];
+