1 #include "simulation.hpp"
6 return (float)rand()/(float)(RAND_MAX);
9 // Helper to account for bounds
10 inline void set(float x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k, float y) {
11 if (i < 0 || i >= CLOUD_DIM_X ||
12 j < 0 || j >= CLOUD_DIM_Y ||
13 k < 0 || k >= CLOUD_DIM_Z)
18 void initClouds(Clouds *cs) {
19 for (int i = 0; i < CLOUD_DIM_X; i++) {
20 for (int j = 0; j < CLOUD_DIM_Y; j++) {
21 for (int k = 0; k < CLOUD_DIM_Z; k++) {
22 cs->act[i][j][k] = randf() < 0.005;
23 cs->cld[i][j][k] = false;
24 cs->hum[i][j][k] = randf() < 0.01;
25 cs->p_ext[i][j][k] = 0.f;
26 cs->p_hum[i][j][k] = 0.f;
27 cs->p_act[i][j][k] = 0.f;
32 for (int k = 0; k < CLOUD_DIM_Y; k++)
33 cs->vy[k] = floor(randf() * 2);
35 // generate ellipsoids of probability
36 const int numEllipsoids = CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z * 0.001;
37 for (int n = 0; n < numEllipsoids; n++) {
38 const float maxSize = 8, minSize = 4;
39 float delta = maxSize - minSize;
40 int width = minSize + randf() * delta, height = minSize + randf() * delta, depth = minSize + randf() * delta;
41 int x = randf() * CLOUD_DIM_X, y = randf() * CLOUD_DIM_Y, z = randf() * CLOUD_DIM_Z;
43 const float maxDist = glm::distance(glm::vec3(0), glm::vec3(0.5,0.5,0.5));
44 for (int i = x; i < x + width; i++) {
45 for (int j = y; j < y + height; j++) {
46 for (int k = z; k < z + depth; k++) {
47 // between [0,0,0] and [1,1,1]
48 glm::vec3 uniformPos = glm::vec3(i - x, j - y, k - z) / glm::vec3(width, height, depth);
50 float dist = glm::distance(uniformPos, glm::vec3(0.5,0.5,0.5)) / maxDist;
51 set(cs->p_ext, i, j, k, P_EXT * dist);
52 set(cs->p_hum, i, j, k, P_HUM * (1.f - dist));
53 set(cs->p_act, i, j, k, P_ACT * (1.f - dist));
60 // Helper to account for bounds
61 inline bool get(bool x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k) {
62 if (i < 0 || i >= CLOUD_DIM_X ||
63 j < 0 || j >= CLOUD_DIM_Y ||
64 k < 0 || k >= CLOUD_DIM_Z)
69 inline bool f_act(Clouds *cs, int i, int j, int k) {
70 return get(cs->act, i + 1, j, k) || get(cs->act, i, j + 1, k)
71 || get(cs->act, i, j, k + 1) || get(cs->act, i - 1, j, k) || get(cs->act, i, j - 1, k)
72 || get(cs->act, i , j, k - 1) || get(cs->act, i - 2, j, k) || get(cs->act, i + 2, j , k)
73 || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2);
76 // scratch for updates on the heap
79 void growth(Clouds *cs) {
82 for (int i = 0; i < CLOUD_DIM_X; i++) {
83 for (int j = 0; j < CLOUD_DIM_Y; j++) {
84 for (int k = 0; k < CLOUD_DIM_Z; k++) {
85 ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
86 ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
87 ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
95 void extinction(Clouds *cs) {
97 for (int i = 0; i < CLOUD_DIM_X; i++) {
98 for (int j = 0; j < CLOUD_DIM_Y; j++) {
99 for (int k = 0; k < CLOUD_DIM_Z; k++) {
100 ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
101 ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
102 ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
109 void advection(Clouds *cs) {
112 for (int i = 0; i < CLOUD_DIM_X; i++) {
113 for (int j = 0; j < CLOUD_DIM_Y; j++) {
114 for (int k = 0; k < CLOUD_DIM_Z; k++) {
116 ncs.hum[i][j][k] = i - v > 0 ? cs->hum[i - v][j][k] : 0;
117 ncs.cld[i][j][k] = i - v > 0 ? cs->cld[i - v][j][k] : 0;
118 ncs.act[i][j][k] = i - v > 0 ? cs->act[i - v][j][k] : 0;
126 /** Weighting function */
127 // TODO: fill this out
128 float w(int ip, int jp, int kp) { return 1; }
130 void calcContDist(Clouds *cls) {
131 const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
132 const float divisor =
133 1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
134 for (int i = 0; i < CLOUD_DIM_X; i++) {
135 for (int j = 0; j < CLOUD_DIM_Y; j++) {
136 for (int k = 0; k < CLOUD_DIM_Z; k++) {
140 for (int tp = -t0; tp <= t0; tp++) {
141 for (int ip = -i0; ip <= i0; ip++) {
142 for (int jp = -j0; jp <= j0; jp++) {
143 for (int kp = -k0; kp <= k0; kp++) {
144 if (i + ip < 0 || i + ip >= CLOUD_DIM_X ||
145 j + jp < 0 || j + jp >= CLOUD_DIM_Y ||
146 k + kp < 0 || k + kp >= CLOUD_DIM_Z)
149 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
155 cls->q[i][j][k] = sum * divisor;
156 if (cls->q[i][j][k] > 1) abort();
162 void stepClouds(Clouds *cs) {