+void advection(Clouds *cs) {
+ Clouds ncs = *cs;
+
+ for (int i = 0; i < CLOUD_DIM_X; i++) {
+ for (int j = 0; j < CLOUD_DIM_Y; j++) {
+ for (int k = 0; k < CLOUD_DIM_Z; k++) {
+ int v = cs->vz[k];
+ ncs.hum[i][j][k] = i - v > 0 ? cs->hum[i - v][j][k] : 0;
+ ncs.cld[i][j][k] = i - v > 0 ? cs->cld[i - v][j][k] : 0;
+ ncs.act[i][j][k] = i - v > 0 ? cs->act[i - v][j][k] : 0;
+ }
+ }
+ }
+
+ *cs = ncs;
+}
+
+/** Weighting function */
+// TODO: fill this out
+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_X; i++) {
+ for (int j = 0; j < CLOUD_DIM_Y; j++) {
+ for (int k = 0; k < CLOUD_DIM_Z; 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_X ||
+ j + jp < 0 || j + jp >= CLOUD_DIM_Y ||
+ k + kp < 0 || k + kp >= CLOUD_DIM_Z)
+ 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();
+ }
+ }
+ }
+}
+