simulation.cpp

   1 #include "simulation.hpp"
   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 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;
  28       }
  29     }
  30   }
  31
  32   for (int k = 0; k < CLOUD_DIM_Y; k++)
  33     cs->vy[k] = floor(randf() * 2);
  34
  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;
  42
  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);
  49
  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));
  54         }
  55       }
  56     }
  57   }
  58 }
  59
  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)
  65     return false;
  66   return x[i][j][k];
  67 }
  68
  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);
  74 }
  75
  76 // scratch for updates on the heap
  77 Clouds ncs;
  78
  79 void growth(Clouds *cs) {
  80   ncs = *cs;
  81
  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);
  88       }
  89     }
  90   }
  91
  92   *cs = ncs;
  93 }
  94
  95 void extinction(Clouds *cs) {
  96   ncs = *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]);
 103       }
 104     }
 105   }
 106   *cs = ncs;
 107 }
 108
 109 void advection(Clouds *cs) {
 110   ncs = *cs;
 111
 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++) {
 115         int v = cs->vy[j];
 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;
 119       }
 120     }
 121   }
 122
 123   *cs = ncs;
 124 }
 125
 126 /** Weighting function */
 127 // TODO: fill this out
 128 float w(int ip, int jp, int kp) { return 1; }
 129
 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++) {
 137         float sum = 0;
 138
 139         // sum
 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)
 147                   continue;
 148
 149                 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
 150               }
 151             }
 152           }
 153         }
 154
 155         cls->q[i][j][k] = sum * divisor;
 156         if (cls->q[i][j][k] > 1) abort();
 157       }
 158     }
 159   }
 160 }
 161
 162 void stepClouds(Clouds *cs) {
 163   growth(cs);
 164   extinction(cs);
 165   advection(cs);
 166   calcContDist(cs);
 167 }