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