X-Git-Url: https://git.lukelau.me/?p=clouds.git;a=blobdiff_plain;f=simulation.cpp;h=2268650cc4250f13be0386a6e0b1ca33b264aa81;hp=360fb07481a2702188c6afd63ea9d34b72c27846;hb=HEAD;hpb=92594780db1f9bb51b97e68457c526fa4287a29e diff --git a/simulation.cpp b/simulation.cpp index 360fb07..2268650 100644 --- a/simulation.cpp +++ b/simulation.cpp @@ -1,4 +1,4 @@ -#include "simulation.h" +#include "simulation.hpp" #include #include @@ -7,23 +7,19 @@ inline float randf() { } // 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] = randf() < 0.01; + 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] = randf() < 0.01; cs->p_ext[i][j][k] = 0.f; @@ -33,28 +29,25 @@ 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); - } - } - } + 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 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, y = randf() * CLOUD_DIM, z = randf() * CLOUD_DIM; - glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2); + 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++) { - float dist = glm::distance(glm::vec3(i,j,k), 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)); @@ -65,10 +58,10 @@ void initClouds(Clouds *cs) { } // 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]; } @@ -80,12 +73,15 @@ inline bool f_act(Clouds *cs, int i, int j, int k) { || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2); } +// scratch for updates on the heap +Clouds ncs; + void growth(Clouds *cs) { - Clouds ncs = *cs; + 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++) { + 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); @@ -97,10 +93,10 @@ void growth(Clouds *cs) { } void extinction(Clouds *cs) { - 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++) { + 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]); @@ -110,16 +106,34 @@ void extinction(Clouds *cs) { *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; i++) { - for (int j = 0; j < CLOUD_DIM; j++) { - for (int k = 0; k < CLOUD_DIM; k++) { + 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 @@ -127,9 +141,9 @@ void calcContDist(Clouds *cls) { 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) + 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]; @@ -148,5 +162,6 @@ void calcContDist(Clouds *cls) { void stepClouds(Clouds *cs) { growth(cs); extinction(cs); + advection(cs); calcContDist(cs); }