Fix ellipsoids in simulation

[clouds.git] / simulation.cpp

diff --git a/simulation.cpp b/simulation.cpp
index 0471f4a20b4cba1dba112237513f940df07051fe..360fb07481a2702188c6afd63ea9d34b72c27846 100644 (file)
--- a/simulation.cpp
+++ b/simulation.cpp
@@ -23,9 +23,9 @@ 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->act[i][j][k] = randf() < 0.01;
         cs->cld[i][j][k] = false;
-        cs->hum[i][j][k] = rand() % 9 == 0;
+        cs->hum[i][j][k] = randf() < 0.01;
         cs->p_ext[i][j][k] = 0.f;
         cs->p_hum[i][j][k] = 0.f;
         cs->p_act[i][j][k] = 0.f;
@@ -33,19 +33,31 @@ 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++) {
+        assert(cs->p_act[i][j][k] == 0.f);
+        assert(cs->p_ext[i][j][k] == 0.f);
+        assert(cs->p_hum[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;
+    const float maxSize = 8, minSize = 4;
+    float delta = maxSize - minSize;
+    int width = minSize + randf() * delta, height = minSize + randf() * delta, depth = minSize + randf() * delta;
     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));
+          float dist = glm::distance(glm::vec3(i,j,k), center) / maxSize;
+          set(cs->p_ext, i, j, k, P_EXT * dist);
+          set(cs->p_hum, i, j, k, P_HUM * (1.f - dist));
+          set(cs->p_act, i, j, k, P_ACT * (1.f - dist));
         }
       }
     }
@@ -69,7 +81,7 @@ inline bool f_act(Clouds *cs, int i, int j, int k) {
 }
 
 void growth(Clouds *cs) {
-  Clouds ncs;
+  Clouds ncs = *cs;
 
   for (int i = 0; i < CLOUD_DIM; i++) {
     for (int j = 0; j < CLOUD_DIM; j++) {
@@ -85,7 +97,7 @@ void growth(Clouds *cs) {
 }
 
 void extinction(Clouds *cs) {
-  Clouds ncs;
+  Clouds ncs = *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++) {
@@ -98,7 +110,43 @@ void extinction(Clouds *cs) {
   *cs = ncs;
 }
 
+/** Weighting function */
+float w(int ip, int jp, int kp) { return 1; }
+
+void calcContDist(Clouds *cls) {
+  const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
+  const float divisor =
+      1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
+  for (int i = 0; i < CLOUD_DIM; i++) {
+    for (int j = 0; j < CLOUD_DIM; j++) {
+      for (int k = 0; k < CLOUD_DIM; k++) {
+        float sum = 0;
+
+        // sum
+        for (int tp = -t0; tp <= t0; tp++) {
+          for (int ip = -i0; ip <= i0; ip++) {
+            for (int jp = -j0; jp <= j0; jp++) {
+              for (int kp = -k0; kp <= k0; kp++) {
+                if (i + ip < 0 || i + ip >= CLOUD_DIM ||
+                    j + jp < 0 || j + jp >= CLOUD_DIM ||
+                    k + kp < 0 || k + kp >= CLOUD_DIM)
+                  continue;
+
+                sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
+              }
+            }
+          }
+        }
+
+        cls->q[i][j][k] = sum * divisor;
+        if (cls->q[i][j][k] > 1) abort();
+      }
+    }
+  }
+}
+
 void stepClouds(Clouds *cs) {
   growth(cs);
   extinction(cs);
+  calcContDist(cs);
 }