-#include "simulation.h"
+#include "simulation.hpp"
#include <cstdlib>
#include <glm/glm.hpp>
}
// 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)
+inline void set(float x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k, float y) {
+ if (i < 0 || i >= CLOUD_DIM_X ||
+ j < 0 || j >= CLOUD_DIM_Y ||
+ k < 0 || k >= CLOUD_DIM_Z)
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;
+ 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++) {
+ cs->act[i][j][k] = randf() < 0.005;
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;
}
}
+ for (int k = 0; k < CLOUD_DIM_Y; k++)
+ cs->vy[k] = floor(randf() * 2);
+
// 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);
+ const int numEllipsoids = CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z * 0.001;
+ for (int n = 0; n < numEllipsoids; n++) {
+ 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_X, y = randf() * CLOUD_DIM_Y, z = randf() * CLOUD_DIM_Z;
+ const float maxDist = glm::distance(glm::vec3(0), glm::vec3(0.5,0.5,0.5));
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));
+ // between [0,0,0] and [1,1,1]
+ glm::vec3 uniformPos = glm::vec3(i - x, j - y, k - z) / glm::vec3(width, height, depth);
+
+ float dist = glm::distance(uniformPos, glm::vec3(0.5,0.5,0.5)) / maxDist;
+ 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));
}
}
}
}
// 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)
+inline bool get(bool x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k) {
+ if (i < 0 || i >= CLOUD_DIM_X ||
+ j < 0 || j >= CLOUD_DIM_Y ||
+ k < 0 || k >= CLOUD_DIM_Z)
return false;
return x[i][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) {
+// scratch for updates on the heap
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++) {
+void growth(Clouds *cs) {
+ 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++) {
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);
}
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 = *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++) {
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 advection(Clouds *cs) {
+ 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->vy[j];
+ 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();
+ }
+ }
+ }
+}
+
void stepClouds(Clouds *cs) {
growth(cs);
extinction(cs);
+ advection(cs);
+ calcContDist(cs);
}
projects / clouds.git / blobdiff