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_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)
  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_X; i++) {
  24     for (int j = 0; j < CLOUD_DIM_Y; j++) {
  25       for (int k = 0; k < CLOUD_DIM_Z; 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.1;
  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_X; i++) {
  37     for (int j = 0; j < CLOUD_DIM_Y; j++) {
  38       for (int k = 0; k < CLOUD_DIM_Z; 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   for (int k = 0; k < CLOUD_DIM_Z; k++)
  46     cs->vz[k] = floor(randf() * 3);
  47
  48   // generate ellipsoids of probability
  49   const int numEllipsoids = CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z * 0.002;
  50   for (int n = 0; n < numEllipsoids; n++) {
  51     const float maxSize = 8, minSize = 4;
  52     float delta = maxSize - minSize;
  53     int width = minSize + randf() * delta, height = minSize + randf() * delta, depth = minSize + randf() * delta;
  54     int x = randf() * CLOUD_DIM_X, y = randf() * CLOUD_DIM_Y, z = randf() * CLOUD_DIM_Z;
  55     glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2);
  56
  57     for (int i = x; i < x + width; i++) {
  58       for (int j = y; j < y + height; j++) {
  59         for (int k = z; k < z + depth; k++) {
  60           float dist = glm::distance(glm::vec3(i,j,k), center) / maxSize;
  61           set(cs->p_ext, i, j, k, P_EXT * dist);
  62           set(cs->p_hum, i, j, k, P_HUM * (1.f - dist));
  63           set(cs->p_act, i, j, k, P_ACT * (1.f - dist));
  64         }
  65       }
  66     }
  67   }
  68 }
  69
  70 // Helper to account for bounds
  71 inline bool get(bool x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k) {
  72   if (i < 0 || i >= CLOUD_DIM_X ||
  73       j < 0 || j >= CLOUD_DIM_Y ||
  74       k < 0 || k >= CLOUD_DIM_Z)
  75     return false;
  76   return x[i][j][k];
  77 }
  78
  79 inline bool f_act(Clouds *cs, int i, int j, int k) {
  80   return get(cs->act, i + 1, j, k) || get(cs->act, i, j + 1, k)
  81     || get(cs->act, i, j, k + 1) || get(cs->act, i - 1, j, k) || get(cs->act, i, j - 1, k)
  82     || get(cs->act, i , j, k - 1) || get(cs->act, i - 2, j, k) || get(cs->act, i + 2, j , k)
  83     || get(cs->act, i, j - 2, k) || get(cs->act, i , j + 2, k) || get(cs->act, i, j, k - 2);
  84 }
  85
  86 void growth(Clouds *cs) {
  87   Clouds ncs = *cs;
  88
  89   for (int i = 0; i < CLOUD_DIM_X; i++) {
  90     for (int j = 0; j < CLOUD_DIM_Y; j++) {
  91       for (int k = 0; k < CLOUD_DIM_Z; k++) {
  92         ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
  93         ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
  94         ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
  95       }
  96     }
  97   }
  98
  99   *cs = ncs;
 100 }
 101
 102 void extinction(Clouds *cs) {
 103   Clouds ncs = *cs;
 104   for (int i = 0; i < CLOUD_DIM_X; i++) {
 105     for (int j = 0; j < CLOUD_DIM_Y; j++) {
 106       for (int k = 0; k < CLOUD_DIM_Z; k++) {
 107         ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
 108         ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
 109         ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
 110       }
 111     }
 112   }
 113   *cs = ncs;
 114 }
 115
 116 void advection(Clouds *cs) {
 117   Clouds ncs = *cs;
 118
 119   for (int i = 0; i < CLOUD_DIM_X; i++) {
 120     for (int j = 0; j < CLOUD_DIM_Y; j++) {
 121       for (int k = 0; k < CLOUD_DIM_Z; k++) {
 122         int v = cs->vz[k];
 123         ncs.hum[i][j][k] = i - v > 0 ? cs->hum[i - v][j][k] : 0;
 124         ncs.cld[i][j][k] = i - v > 0 ? cs->cld[i - v][j][k] : 0;
 125         ncs.act[i][j][k] = i - v > 0 ? cs->act[i - v][j][k] : 0;
 126       }
 127     }
 128   }
 129
 130   *cs = ncs;
 131 }
 132
 133 /** Weighting function */
 134 // TODO: fill this out
 135 float w(int ip, int jp, int kp) { return 1; }
 136
 137 void calcContDist(Clouds *cls) {
 138   const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
 139   const float divisor =
 140       1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
 141   for (int i = 0; i < CLOUD_DIM_X; i++) {
 142     for (int j = 0; j < CLOUD_DIM_Y; j++) {
 143       for (int k = 0; k < CLOUD_DIM_Z; k++) {
 144         float sum = 0;
 145
 146         // sum
 147         for (int tp = -t0; tp <= t0; tp++) {
 148           for (int ip = -i0; ip <= i0; ip++) {
 149             for (int jp = -j0; jp <= j0; jp++) {
 150               for (int kp = -k0; kp <= k0; kp++) {
 151                 if (i + ip < 0 || i + ip >= CLOUD_DIM_X ||
 152                     j + jp < 0 || j + jp >= CLOUD_DIM_Y ||
 153                     k + kp < 0 || k + kp >= CLOUD_DIM_Z)
 154                   continue;
 155
 156                 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
 157               }
 158             }
 159           }
 160         }
 161
 162         cls->q[i][j][k] = sum * divisor;
 163         if (cls->q[i][j][k] > 1) abort();
 164       }
 165     }
 166   }
 167 }
 168
 169 void stepClouds(Clouds *cs) {
 170   growth(cs);
 171   extinction(cs);
 172   advection(cs);
 173   calcContDist(cs);
 174 }