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] = 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;
  32       }
  33     }
  34   }
  35
  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);
  42
  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));
  49         }
  50       }
  51     }
  52   }
  53 }
  54
  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)
  60     return false;
  61   return x[i][j][k];
  62 }
  63
  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);
  69 }
  70
  71 void growth(Clouds *cs) {
  72   Clouds ncs;
  73
  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);
  80       }
  81     }
  82   }
  83
  84   *cs = ncs;
  85 }
  86
  87 void extinction(Clouds *cs) {
  88   Clouds ncs;
  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]);
  95       }
  96     }
  97   }
  98   *cs = ncs;
  99 }
 100
 101 void stepClouds(Clouds *cs) {
 102   growth(cs);
 103   extinction(cs);
 104 }