Start simulation

[clouds.git] / simulation.cpp

diff --git a/simulation.cpp b/simulation.cpp
new file mode 100644 (file)
index 0000000..0471f4a
--- /dev/null
+++ b/simulation.cpp
@@ -0,0 +1,104 @@
+#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);
+}