simulation.cpp

   1 #include "simulation.h"
   2 #include <cstdlib>
   3 #include <glm/glm.hpp>
   4
   5 inline float randf() {
   6   return  (float)rand()/(float)(RAND_MAX);
   7 }
   8
   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)
  14     return;
  15   x[i][j][k] = y;
  16 }
  17
  18 #define P_EXT 0.1
  19 #define P_HUM 0.1
  20 #define P_ACT 0.001
  21
  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] = randf() < 0.01;
  27         cs->cld[i][j][k] = false;
  28         cs->hum[i][j][k] = randf() < 0.01;
  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;
  32       }
  33     }
  34   }
  35
  36   for (int i = 0; i < CLOUD_DIM; i++) {
  37     for (int j = 0; j < CLOUD_DIM; j++) {
  38       for (int k = 0; k < CLOUD_DIM; k++) {
  39         assert(cs->p_act[i][j][k] == 0.f);
  40         assert(cs->p_ext[i][j][k] == 0.f);
  41         assert(cs->p_hum[i][j][k] == 0.f);
  42       }
  43     }
  44   }
  45
  46   // generate ellipsoids of probability
  47   for (int n = 0; n < 6; n++) {
  48     const float maxSize = 8, minSize = 4;
  49     float delta = maxSize - minSize;
  50     int width = minSize + randf() * delta, height = minSize + randf() * delta, depth = minSize + randf() * delta;
  51     int x = randf() * CLOUD_DIM, y = randf() * CLOUD_DIM, z = randf() * CLOUD_DIM;
  52     glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2);
  53
  54     for (int i = x; i < x + width; i++) {
  55       for (int j = y; j < y + height; j++) {
  56         for (int k = z; k < z + depth; k++) {
  57           float dist = glm::distance(glm::vec3(i,j,k), center) / maxSize;
  58           set(cs->p_ext, i, j, k, P_EXT * dist);
  59           set(cs->p_hum, i, j, k, P_HUM * (1.f - dist));
  60           set(cs->p_act, i, j, k, P_ACT * (1.f - dist));
  61         }
  62       }
  63     }
  64   }
  65 }
  66
  67 // Helper to account for bounds
  68 inline bool get(bool x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k) {
  69   if (i < 0 || i >= CLOUD_DIM ||
  70       j < 0 || j >= CLOUD_DIM ||
  71       k < 0 || k >= CLOUD_DIM)
  72     return false;
  73   return x[i][j][k];
  74 }
  75
  76 inline bool f_act(Clouds *cs, int i, int j, int k) {
  77   return get(cs->act, i + 1, j, k) || get(cs->act, i, j + 1, k)
  78     || get(cs->act, i, j, k + 1) || get(cs->act, i - 1, j, k) || get(cs->act, i, j - 1, k)
  79     || get(cs->act, i , j, k - 1) || get(cs->act, i - 2, j, k) || get(cs->act, i + 2, j , k)
  80     || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2);
  81 }
  82
  83 void growth(Clouds *cs) {
  84   Clouds ncs = *cs;
  85
  86   for (int i = 0; i < CLOUD_DIM; i++) {
  87     for (int j = 0; j < CLOUD_DIM; j++) {
  88       for (int k = 0; k < CLOUD_DIM; k++) {
  89         ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
  90         ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
  91         ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
  92       }
  93     }
  94   }
  95
  96   *cs = ncs;
  97 }
  98
  99 void extinction(Clouds *cs) {
 100   Clouds ncs = *cs;
 101   for (int i = 0; i < CLOUD_DIM; i++) {
 102     for (int j = 0; j < CLOUD_DIM; j++) {
 103       for (int k = 0; k < CLOUD_DIM; k++) {
 104         ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
 105         ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
 106         ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
 107       }
 108     }
 109   }
 110   *cs = ncs;
 111 }
 112
 113 /** Weighting function */
 114 float w(int ip, int jp, int kp) { return 1; }
 115
 116 void calcContDist(Clouds *cls) {
 117   const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
 118   const float divisor =
 119       1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
 120   for (int i = 0; i < CLOUD_DIM; i++) {
 121     for (int j = 0; j < CLOUD_DIM; j++) {
 122       for (int k = 0; k < CLOUD_DIM; k++) {
 123         float sum = 0;
 124
 125         // sum
 126         for (int tp = -t0; tp <= t0; tp++) {
 127           for (int ip = -i0; ip <= i0; ip++) {
 128             for (int jp = -j0; jp <= j0; jp++) {
 129               for (int kp = -k0; kp <= k0; kp++) {
 130                 if (i + ip < 0 || i + ip >= CLOUD_DIM ||
 131                     j + jp < 0 || j + jp >= CLOUD_DIM ||
 132                     k + kp < 0 || k + kp >= CLOUD_DIM)
 133                   continue;
 134
 135                 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
 136               }
 137             }
 138           }
 139         }
 140
 141         cls->q[i][j][k] = sum * divisor;
 142         if (cls->q[i][j][k] > 1) abort();
 143       }
 144     }
 145   }
 146 }
 147
 148 void stepClouds(Clouds *cs) {
 149   growth(cs);
 150   extinction(cs);
 151   calcContDist(cs);
 152 }