Continuous distribution

[clouds.git] / simulation.cpp

diff --git a/simulation.cpp b/simulation.cpp
index 0471f4a20b4cba1dba112237513f940df07051fe..e022d12b1b921901b42dcceaa50466a0e96b439f 100644 (file)
--- a/simulation.cpp
+++ b/simulation.cpp
@@ -98,7 +98,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);
 }