1 #include "simulation.h"
6 return (float)rand()/(float)(RAND_MAX);
9 // Helper to account for bounds
10 inline void set(float x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k, float y) {
11 if (i < 0 || i >= CLOUD_DIM ||
12 j < 0 || j >= CLOUD_DIM ||
13 k < 0 || k >= CLOUD_DIM)
22 void initClouds(Clouds *cs) {
23 for (int i = 0; i < CLOUD_DIM; i++) {
24 for (int j = 0; j < CLOUD_DIM; j++) {
25 for (int k = 0; k < CLOUD_DIM; k++) {
26 cs->act[i][j][k] = rand() % 8 == 0;
27 cs->cld[i][j][k] = false;
28 cs->hum[i][j][k] = rand() % 9 == 0;
29 cs->p_ext[i][j][k] = 0.f;
30 cs->p_hum[i][j][k] = 0.f;
31 cs->p_act[i][j][k] = 0.f;
36 // generate ellipsoids of probability
37 for (int n = 0; n < 6; n++) {
38 const float maxSize = 5;
39 int width = randf() * maxSize, height = randf() * maxSize, depth = randf() * maxSize;
40 int x = randf() * CLOUD_DIM, y = randf() * CLOUD_DIM, z = randf() * CLOUD_DIM;
41 glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2);
43 for (int i = x; i < x + width; i++) {
44 for (int j = y; j < y + height; j++) {
45 for (int k = z; k < z + depth; k++) {
46 set(cs->p_ext, x, y, z, P_EXT * glm::distance(glm::vec3(x,y,z), center) / maxSize);
47 set(cs->p_hum, x, y, z, P_HUM * (1.f - glm::distance(glm::vec3(x,y,z), center) / maxSize));
48 set(cs->p_act, x, y, z, P_ACT * (1.f - glm::distance(glm::vec3(x,y,z), center) / maxSize));
55 // Helper to account for bounds
56 inline bool get(bool x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k) {
57 if (i < 0 || i >= CLOUD_DIM ||
58 j < 0 || j >= CLOUD_DIM ||
59 k < 0 || k >= CLOUD_DIM)
64 inline bool f_act(Clouds *cs, int i, int j, int k) {
65 return get(cs->act, i + 1, j, k) || get(cs->act, i, j + 1, k)
66 || get(cs->act, i, j, k + 1) || get(cs->act, i - 1, j, k) || get(cs->act, i, j - 1, k)
67 || get(cs->act, i , j, k - 1) || get(cs->act, i - 2, j, k) || get(cs->act, i + 2, j , k)
68 || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2);
71 void growth(Clouds *cs) {
74 for (int i = 0; i < CLOUD_DIM; i++) {
75 for (int j = 0; j < CLOUD_DIM; j++) {
76 for (int k = 0; k < CLOUD_DIM; k++) {
77 ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
78 ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
79 ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
87 void extinction(Clouds *cs) {
89 for (int i = 0; i < CLOUD_DIM; i++) {
90 for (int j = 0; j < CLOUD_DIM; j++) {
91 for (int k = 0; k < CLOUD_DIM; k++) {
92 ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
93 ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
94 ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
101 /** Weighting function */
102 float w(int ip, int jp, int kp) { return 1; }
104 void calcContDist(Clouds *cls) {
105 const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
106 const float divisor =
107 1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
108 for (int i = 0; i < CLOUD_DIM; i++) {
109 for (int j = 0; j < CLOUD_DIM; j++) {
110 for (int k = 0; k < CLOUD_DIM; k++) {
114 for (int tp = -t0; tp <= t0; tp++) {
115 for (int ip = -i0; ip <= i0; ip++) {
116 for (int jp = -j0; jp <= j0; jp++) {
117 for (int kp = -k0; kp <= k0; kp++) {
118 if (i + ip < 0 || i + ip >= CLOUD_DIM ||
119 j + jp < 0 || j + jp >= CLOUD_DIM ||
120 k + kp < 0 || k + kp >= CLOUD_DIM)
123 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
129 cls->q[i][j][k] = sum * divisor;
130 if (cls->q[i][j][k] > 1) abort();
136 void stepClouds(Clouds *cs) {