+#include "simulation.h"
+#include <cstdlib>
+#include <glm/glm.hpp>
+
+inline float randf() {
+ return (float)rand()/(float)(RAND_MAX);
+}
+
+// Helper to account for bounds
+inline void set(float x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k, float y) {
+ if (i < 0 || i >= CLOUD_DIM ||
+ j < 0 || j >= CLOUD_DIM ||
+ k < 0 || k >= CLOUD_DIM)
+ return;
+ x[i][j][k] = y;
+}
+
+#define P_EXT 0.1
+#define P_HUM 0.1
+#define P_ACT 0.001
+
+void initClouds(Clouds *cs) {
+ for (int i = 0; i < CLOUD_DIM; i++) {
+ for (int j = 0; j < CLOUD_DIM; j++) {
+ for (int k = 0; k < CLOUD_DIM; k++) {
+ cs->act[i][j][k] = rand() % 8 == 0;
+ cs->cld[i][j][k] = false;
+ cs->hum[i][j][k] = rand() % 9 == 0;
+ cs->p_ext[i][j][k] = 0.f;
+ cs->p_hum[i][j][k] = 0.f;
+ cs->p_act[i][j][k] = 0.f;
+ }
+ }
+ }
+
+ // generate ellipsoids of probability
+ for (int n = 0; n < 6; n++) {
+ const float maxSize = 5;
+ int width = randf() * maxSize, height = randf() * maxSize, depth = randf() * maxSize;
+ int x = randf() * CLOUD_DIM, y = randf() * CLOUD_DIM, z = randf() * CLOUD_DIM;
+ glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2);
+
+ for (int i = x; i < x + width; i++) {
+ for (int j = y; j < y + height; j++) {
+ for (int k = z; k < z + depth; k++) {
+ set(cs->p_ext, x, y, z, P_EXT * glm::distance(glm::vec3(x,y,z), center) / maxSize);
+ set(cs->p_hum, x, y, z, P_HUM * (1.f - glm::distance(glm::vec3(x,y,z), center) / maxSize));
+ set(cs->p_act, x, y, z, P_ACT * (1.f - glm::distance(glm::vec3(x,y,z), center) / maxSize));
+ }
+ }
+ }
+ }
+}
+
+// Helper to account for bounds
+inline bool get(bool x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k) {
+ if (i < 0 || i >= CLOUD_DIM ||
+ j < 0 || j >= CLOUD_DIM ||
+ k < 0 || k >= CLOUD_DIM)
+ return false;
+ return x[i][j][k];
+}
+
+inline bool f_act(Clouds *cs, int i, int j, int k) {
+ return get(cs->act, i + 1, j, k) || get(cs->act, i, j + 1, k)
+ || get(cs->act, i, j, k + 1) || get(cs->act, i - 1, j, k) || get(cs->act, i, j - 1, k)
+ || get(cs->act, i , j, k - 1) || get(cs->act, i - 2, j, k) || get(cs->act, i + 2, j , k)
+ || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2);
+}
+
+void growth(Clouds *cs) {
+ Clouds ncs;
+
+ for (int i = 0; i < CLOUD_DIM; i++) {
+ for (int j = 0; j < CLOUD_DIM; j++) {
+ for (int k = 0; k < CLOUD_DIM; k++) {
+ ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
+ ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
+ ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
+ }
+ }
+ }
+
+ *cs = ncs;
+}
+
+void extinction(Clouds *cs) {
+ Clouds ncs;
+ for (int i = 0; i < CLOUD_DIM; i++) {
+ for (int j = 0; j < CLOUD_DIM; j++) {
+ for (int k = 0; k < CLOUD_DIM; k++) {
+ ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
+ ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
+ ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
+ }
+ }
+ }
+ *cs = ncs;
+}
+
+void stepClouds(Clouds *cs) {
+ growth(cs);
+ extinction(cs);
+}