typedef struct {
  HMap *hmap;
  HMap *wmap;
} Water;



void waterInit(Water *w, HMap *hmap, HMap *wmap) {
  w->hmap = hmap;
  w->wmap = wmap;
  for(int i = 0; i < HMAP_SIZE; ++i)
    for(int j = 0; j < HMAP_SIZE; ++j)
      w->wmap->map[i][j] = w->hmap->map[i][j];
}


void waterRain(Water *w, double k) {
  for(int i = 0; i < HMAP_SIZE; ++i)
    for(int j = 0; j < HMAP_SIZE; ++j)
      w->wmap->map[i][j] += k;
}


void waterFall(Water *w, double wk, double hk) {
  static const int di[4] = { -1, 1,  0, 0 };
  static const int dj[4] = {  0, 0, -1, 1 };

  hmapCopy(w->hmap);
  hmapCopy(w->wmap);
  for(int i = 0; i < HMAP_SIZE; ++i) {
    for(int j = 0; j < HMAP_SIZE; ++j) {
      double dz = w->wmap->old[i][j] - w->hmap->old[i][j];
      dz *= 0.25;
      if (dz > 0) {
        double sum = 0;
        double v[4] = {};
        for(int k = 0; k < 4; ++k) {
          int ii = ((i+di[k]) + HMAP_SIZE)%HMAP_SIZE;
          int jj = ((j+dj[k]) + HMAP_SIZE)%HMAP_SIZE;
          double vv = w->wmap->old[i][j] - w->wmap->old[ii][jj];
          if (vv < 0) vv = 0;
          //vv = vv*vv*vv;
          v[k] = vv;
          sum += vv;
        }
        double s = wk*sum;
        if (s > dz) s = dz;
        if (sum > 1e-10) {
          for(int k = 0; k < 4; ++k) {
            int ii = ((i+di[k]) + HMAP_SIZE)%HMAP_SIZE;
            int jj = ((j+dj[k]) + HMAP_SIZE)%HMAP_SIZE;
            double vv = v[k]/sum*s;
            w->wmap->map[ii][jj] += v[k]/sum*s;
            w->wmap->map[i][j] -= vv;
            double dh = w->hmap->old[i][j] - w->hmap->map[ii][jj];
            dh *= 0.25;
            if (dh > 0) {
              vv = vv*hk;
              if (vv > dh) vv = dh;
              w->hmap->map[ii][jj] += vv;
              w->hmap->map[i][j] -= vv;
            }
          }
        }
      }
    }
  }
}


double waterEvaporation(Water *w, double k) {
  double sum;
  for(int i = 0; i < HMAP_SIZE; ++i) {
    for(int j = 0; j < HMAP_SIZE; ++j) {
      double dz = w->wmap->map[i][j] - w->hmap->map[i][j];
      if (dz < 0) dz = 0;
      if (dz > k) dz = k;
      sum += dz;
      w->wmap->map[i][j] -= dz;
    }
  }
  return sum;
}


void waterCycle(Water *w, double fallK, double erosionK, double evaK) {
  double sum = waterEvaporation(w, evaK);
  waterRain(w, sum/(HMAP_SIZE*HMAP_SIZE));
  for(int i = 0; i < 10; ++i) {
    waterFall(w, fallK/10, erosionK/10);
    printf(".");
    fflush(stdout);
  }
  printf("\n");
}