Water simulation using height fields and flow maps

Water.cpp

#include "Water.h"
#include <GLFW/glfw3.h>

// --------------------------------------------------------------------------------

int water_get_num_faces(const SMikkTSpaceContext* ctx) {
  Water* w = static_cast<Water*>(ctx->m_pUserData);
  return w->faces.size();
}

void water_get_position(const SMikkTSpaceContext* ctx, float pos[], const int iface, const int ivert) {
  Water* w = static_cast<Water*>(ctx->m_pUserData);
  vec3& v = w->vertices[w->faces[iface].indices[ivert]].pos;
  pos[0] = v.x;
  pos[1] = v.y;
  pos[2] = v.z;
}

void water_get_normal(const SMikkTSpaceContext* ctx, float norm[], const int iface, const int ivert) {
  Water* w = static_cast<Water*>(ctx->m_pUserData);
  vec3& n = w->vertices[w->faces[iface].indices[ivert]].norm;
  norm[0] = n.x;
  norm[1] = n.y;
  norm[2] = n.z;
}

void water_get_texcoord(const SMikkTSpaceContext* ctx, float tex[], const int iface, const int ivert) {
  Water* w = static_cast<Water*>(ctx->m_pUserData);
  vec2& t = w->vertices[w->faces[iface].indices[ivert]].tex;
  tex[0] = t.x;
  tex[1] = t.y;
}

void water_set_tspace(const SMikkTSpaceContext* ctx, const float tangent[], const float sign, const int iface, int ivert) {
}

// --------------------------------------------------------------------------------

Water::Water() 
  :vbo(0)
   //  ,vbo_els(0)
  ,vao(0)
  ,vert(0)
  ,frag(0)
  ,prog(0)
  ,bytes_allocated(0)
  ,flow_tex(0)
  ,norm_tex(0)
  ,noise_tex(0)
  ,diffuse_tex(0)
  ,cube_tex(0)
{
  memset(u0, 0x00, sizeof(u0));
  memset(u1, 0x00, sizeof(u1));
  memset(v, 0x00, sizeof(v));
}

bool Water::setup(int w, int h) {
  pm.perspective(65.0f, 1280.0f/720.0f, 0.01, 100.0f);
  //vm.lookAt(vec3(0.0, 25.0f, 0.0), vec3(0.0,0.0,0.0), vec3(0.0f, 1.0f, 0.0f));

  // TEMP
  //vm.lookAt(vec3(0.0, 15.0f, 0.0), vec3(0.0,0.0,0.0), vec3(0.0f, 1.0f, 0.0f));
  vm.translate(0, -0.4, -5);
  vm.print();
  // TEMP

  vert = rx_create_shader(GL_VERTEX_SHADER, WATER_VS);
  frag = rx_create_shader(GL_FRAGMENT_SHADER, WATER_FS);
  prog = rx_create_program(vert, frag);
  glBindAttribLocation(prog, 0, "a_pos");
  glBindAttribLocation(prog, 1, "a_tex");
  glBindAttribLocation(prog, 2, "a_norm");
  glLinkProgram(prog);
  rx_print_shader_link_info(prog);

  glUseProgram(prog);
  glUniform1i(glGetUniformLocation(prog, "u_norm_tex"), 0);
  glUniform1i(glGetUniformLocation(prog, "u_flow_tex"), 1);
  glUniform1i(glGetUniformLocation(prog, "u_noise_tex"), 2);
  glUniform1i(glGetUniformLocation(prog, "u_diffuse_tex"), 3);
  glUniform1i(glGetUniformLocation(prog, "u_cube_tex"), 4);

  glGenVertexArrays(1, &vao);
  glBindVertexArray(vao);
  glGenBuffers(1, &vbo);
  glBindBuffer(GL_ARRAY_BUFFER, vbo);
  glEnableVertexAttribArray(0); // pos
  glEnableVertexAttribArray(1); // tex
  glEnableVertexAttribArray(2); // norm
  glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(WaterVertex), (GLvoid*)0);
  glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, sizeof(WaterVertex), (GLvoid*)12);
  glVertexAttribPointer(2, 3, GL_FLOAT, GL_FALSE, sizeof(WaterVertex), (GLvoid*)20);

  createVertices();
  updateVertices();
  
  //  glGenBuffers(1, &vbo_els);
  //  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbo_els);
  //  glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint) * indices.size(), &indices[0], GL_STATIC_DRAW);

  norm_tex = createTexture("water_normals.png");
  flow_tex = createTexture("water_flow.png");
  noise_tex = createTexture("water_noise.png");
  diffuse_tex = createTexture("water_diffuse.png");
  cube_tex = createCubeMap("debug");

  if(!setupHeightField()) {
    printf("Error: cannot setup the heightfield.\n");
    return false;
  }


  mikk_context.m_pUserData = this;
  mikk_context.m_pInterface = &mikk_interface;
  mikk_interface.m_getNumFaces = water_get_num_faces;
  mikk_interface.m_getNumVerticesOfFace = water_get_num_vertices_of_face;
  mikk_interface.m_getPosition = water_get_position;
  mikk_interface.m_getNormal = water_get_normal;
  mikk_interface.m_getTexCoord = water_get_texcoord;
  mikk_interface.m_setTSpaceBasic = water_set_tspace;
   

  return true;
}

bool Water::setupHeightField() {

  for(int j = 0; j < FIELD_N; ++j) {
    for(int i = 0; i < FIELD_N; ++i) {
      int dx = j * FIELD_N + i;
      if(i > 20 && i < 50 && j > 10 && j < 40) {
        // u0[dx] = sin(float(j)/FIELD_N) * 4.3; // rx_random(0.1f, 0.2f);
        u0[dx]= -2.9f;
      }
      v[dx] = 0.0f;
    }
  }


  return true;
}


void Water::testFlowFieldForce() {
  setupHeightField();
}

// See here for a version which applies a force: 
// https://gist.github.com/roxlu/0cb3f8d2cf17b2bbef66
void Water::addDrop(int col, int row) {
  assert(col < FIELD_N);
  assert(row < FIELD_N);
  setupHeightField();
}

void Water::createVertices() {

  int grid_size = FIELD_N;
  float size = 0.2;
  float hs = size * (grid_size-1) * 0.5;

  tmp_vertices.clear();
  tmp_vertices.assign(grid_size * grid_size, WaterVertex());

  //std::vector<WaterVertex> verts;
  //verts.assign(grid_size * grid_size, WaterVertex());

  for(int j = 0; j < grid_size; ++j) {
    for(int i = 0; i < grid_size; ++i) {
      int dx = j * grid_size + i;
      WaterVertex& v = tmp_vertices[dx];
      v.pos.set(-hs + (i * size), 0.0, -hs + (j * size));
      v.tex.set(float(i)/(grid_size-1), float(j)/(grid_size-1));
      v.norm.set(0.0f, 1.0f, 0.0f);
    }
  }

}

void Water::updateVertices() {

  // triangulate
  vertices.clear();
  faces.clear();
  indices.clear();
  
  for(int j = 0; j < FIELD_N - 1; ++j) {
    for(int i = 0; i < FIELD_N - 1; ++i) {
      GLuint a = ((j + 0) * FIELD_N) + (i + 0);  // bottom left
      GLuint b = ((j + 0) * FIELD_N) + (i + 1);  // bottom right
      GLuint c = ((j + 1) * FIELD_N) + (i + 1);  // top right 
      GLuint d = ((j + 1) * FIELD_N) + (i + 0);  // top left
      
      int tri = vertices.size();
      vertices.push_back(tmp_vertices[a]);
      vertices.push_back(tmp_vertices[b]);
      vertices.push_back(tmp_vertices[c]);
      vertices.push_back(tmp_vertices[a]);
      vertices.push_back(tmp_vertices[c]);
      vertices.push_back(tmp_vertices[d]);

      WaterFace fa(tri,tri+1,tri+2);
      WaterFace fb(tri+3,tri+4,tri+5);
      faces.push_back(fa);
      faces.push_back(fb);
    }
  }


  // update GL 
  glBindBuffer(GL_ARRAY_BUFFER, vbo);

  size_t needed = sizeof(WaterVertex) * vertices.size();
  if(needed > bytes_allocated) {
    glBufferData(GL_ARRAY_BUFFER, needed, vertices[0].pos.ptr(), GL_STATIC_DRAW);
    bytes_allocated = needed;
  }
  else {
    glBufferSubData(GL_ARRAY_BUFFER, 0, needed, vertices[0].pos.ptr()); 
  }
}


void Water::applyHeightFieldToWater() {

  if(!vertices.size() || vertices.size() < FIELD_NN) {
    printf("Error: not enough vetices.\n");
    return;
  }

  for(int j = 0; j < FIELD_N; ++j) {
    for(int i = 0; i < FIELD_N; ++i) {
      int dx = j * FIELD_N + i;
      tmp_vertices[dx].pos.y = u0[dx];
    }
  }
}

void Water::calculateNormals() {
      
  return;
  
  assert(faces.size());

  // @todo - fix borders
  for(int i = 1; i < FIELD_N-1; ++i) {
    for(int j = 1; j < FIELD_N-1; ++j) {
      int dx_current = ((j + 0) * FIELD_N) + (i + 0);
      int dx_left    = ((j + 0) * FIELD_N) + (i - 1);
      int dx_right   = ((j + 0) * FIELD_N) + (i + 1);
      int dx_bottom  = ((j - 1) * FIELD_N) + (i + 0);
      int dx_top     = ((j + 1) * FIELD_N) + (i + 0);

      vec3& n  = vertices[dx_current].pos;
      vec3& a  = vertices[dx_left].pos;
      vec3& b  = vertices[dx_right].pos;
      vec3& c  = vertices[dx_top].pos;
      vec3& d  = vertices[dx_bottom].pos;

      vec3& center  = vertices[dx_current].pos;
      vec3& left  = vertices[dx_left].pos;
      vec3& right  = vertices[dx_right].pos;
      vec3& top  = vertices[dx_top].pos;
      vec3& bottom  = vertices[dx_bottom].pos;
      
      vec3 to_top = top - center;
      vec3 to_right = right - center;
      vec3 to_bottom = bottom - center;
      vec3 to_left = left - center;
      
      vec3 c1 = (cross(to_top, to_right));
      vec3 c2 = (cross(to_left, to_top));
      vec3 c3 = (cross(to_top, to_bottom));
      vec3 c4 = (cross(to_bottom, to_right));
      
      vertices[dx_current].norm = (c1 + c2 + c3 + c4) * 0.25;
    }
  }
#if 0
  // Calculate the normals per face
  for(size_t i = 0; i < faces.size(); ++i) {
    WaterFace& face = faces[i];
    vec3& pos_a = vertices[face.a].pos;
    vec3& pos_b = vertices[face.b].pos;
    vec3& pos_c = vertices[face.c].pos;
    vec3 ab = pos_a - pos_b;
    vec3 ac = pos_c - pos_a;
    face.norm = normalized(cross(ab, ac));
  }

  // Average normals 
  // @todo

  // Set normals 
  // @todo - wrong implementation now
  for(size_t i = 0; i < faces.size(); ++i) {
    WaterFace& face = faces[i];
    vertices[face.a].norm = face.norm;
    vertices[face.b].norm = face.norm;
    vertices[face.c].norm = face.norm;
  }
#endif
}

// takes too long
void Water::calculateTangents() {
  return;
  double n = glfwGetTime();
  genTangSpaceDefault(&mikk_context);
  double dt = glfwGetTime() - n;
  printf("Took: %f\n", dt);
}


void Water::update(float dt) {
  updateHeightField(dt);
  applyHeightFieldToWater();
  calculateNormals();
  calculateTangents();
  updateVertices();
}

void Water::updateHeightField(float dt) {
  float c = 1.0;
  float max_c = (1.0 / dt);
  if(c > max_c) {
    printf("Warning: invalid C value\n");
    return;
  }

  dt = dt * 4.9;
  for(int j = 1; j < FIELD_N - 1; ++j) {
    for(int i = 1; i < FIELD_N - 1; ++i) {
      int current = ((j + 0) * FIELD_N) + (i + 0);
      int right   = ((j + 0) * FIELD_N) + (i + 1);
      int left    = ((j + 0) * FIELD_N) + (i - 1);
      int top     = ((j + 1) * FIELD_N) + (i + 0);
      int bottom  = ((j - 1) * FIELD_N) + (i + 0);
      float f = (c*c) * ((u0[right] + u0[left] + u0[bottom] + u0[top])) - (4.0 * u0[current]) ;
      if(f > 0.1) {
        f = 0.1;
      }
      else if(f < -0.1) {
        f = -0.1;
      }
   
      v[current] = v[current] + f * dt;
      u1[current] = u0[current] + v[current] * dt;
      v[current] *= 0.995;
   
    }
  }
  memcpy(u0, u1, sizeof(u0));
}

void Water::draw(){
  static float t = 0.0;

  float time0 = ((fmod(t, 1.0)));
  float time1 = ((fmod(t + 0.5, 1.0)));
  t += 0.001;

  glEnable(GL_DEPTH_TEST);
  glUseProgram(prog);
  glBindVertexArray(vao);

  glActiveTexture(GL_TEXTURE0);
  glBindTexture(GL_TEXTURE_2D, norm_tex);

  glActiveTexture(GL_TEXTURE1);
  glBindTexture(GL_TEXTURE_2D, flow_tex);

  glActiveTexture(GL_TEXTURE2);
  glBindTexture(GL_TEXTURE_2D, noise_tex);

  glActiveTexture(GL_TEXTURE3);
  glBindTexture(GL_TEXTURE_2D, diffuse_tex);

  glActiveTexture(GL_TEXTURE4);
  glBindTexture(GL_TEXTURE_CUBE_MAP, cube_tex);

  glUniform1f(glGetUniformLocation(prog, "u_time"), t);  
  glUniform1f(glGetUniformLocation(prog, "u_time0"), time0);
  glUniform1f(glGetUniformLocation(prog, "u_time1"), time1);
  glUniformMatrix4fv(glGetUniformLocation(prog, "u_pm"), 1, GL_FALSE, pm.ptr());
  glUniformMatrix4fv(glGetUniformLocation(prog, "u_vm"), 1, GL_FALSE, vm.ptr());
  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  // glDrawElements(GL_TRIANGLES, indices.size(), GL_UNSIGNED_INT, NULL);
  glDrawArrays(GL_TRIANGLES, 0, vertices.size());
}


GLuint Water::createTexture(std::string filename) {

  int w, h, n;
  unsigned char* pix;

  if(!rx_load_png(rx_to_data_path(filename), &pix, w, h, n)) {
    printf("Error: cannot find: %s\n", filename.c_str());
    return 0;
  }

  GLuint tex;
  GLenum format = GL_RGB;

  if(n == 4) {
    format = GL_RGBA;
  }

  glGenTextures(1, &tex);
  glBindTexture(GL_TEXTURE_2D, tex);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, w, h, 0, format, GL_UNSIGNED_BYTE, pix);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  printf("Texture, w: %d, h: %d, n: %d\n", w, h, n);

  delete[] pix;
  pix = NULL;
  return tex;
}


GLuint Water::createCubeMap(std::string prefix) {

  std::string names[] = {"_xpos.png", "_xneg.png", "_ypos.png", "_yneg.png", "_zpos.png", "_zneg.png" };

  GLenum cube[] = {
    GL_TEXTURE_CUBE_MAP_POSITIVE_X,
    GL_TEXTURE_CUBE_MAP_NEGATIVE_X,
    GL_TEXTURE_CUBE_MAP_POSITIVE_Y,
    GL_TEXTURE_CUBE_MAP_NEGATIVE_Y,
    GL_TEXTURE_CUBE_MAP_POSITIVE_Z,
    GL_TEXTURE_CUBE_MAP_NEGATIVE_Z
  };

  GLuint tex = 0;
  glGenTextures(1, &tex);
  glBindTexture(GL_TEXTURE_CUBE_MAP, tex);
  glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);


  int w, h, n = 0;
  unsigned char* pix = NULL;
  GLenum format = GL_RGBA;

  for(int i = 0; i < 6; ++i) {
    std::string filepath = rx_to_data_path(prefix +names[i]);
    
    if(!rx_load_png(filepath, &pix, w, h, n)) {
      printf("Error: cannot load: %s\n", filepath.c_str());
      ::exit(EXIT_FAILURE);
    }

    if(n == 3) {
      format = GL_RGB;
    }
    else if(n == 4) {
      format = GL_RGBA;
    }

    printf("%s, %d x %d x %d\n", filepath.c_str(), w, h, n);

    glTexImage2D(cube[i], 0, GL_RGBA, w, h, 0, format, GL_UNSIGNED_BYTE, pix);

    delete[] pix;
    pix = NULL;
    
  }
  //  glGenerateMipmap(GL_TEXTURE_CUBE_MAP);

  return tex;
}

Water.h

#ifndef WATER_H
#define WATER_H

// Guild Wars                           - http://www.youtube.com/watch?v=GlwvxNuhCIs
// Guild Wars                           - http://www.youtube.com/watch?v=8tRlHa4E32Y
// Guild Wars NICE!                     - http://www.youtube.com/watch?feature=player_detailpage&v=Xm4h2l2CkwU#t=24
// Real time flow map creation          - http://vimeo.com/57182020
// From Dust, nice shores + water       - http://www.youtube.com/watch?v=rTrPyMYdjJM
// From Dust,                           - http://www.youtube.com/watch?v=LVjkgPlKu1g
// From Dust (same as above, water drop - http://www.youtube.com/watch?feature=player_detailpage&v=LVjkgPlKu1g#t=592


// OpenGL, water color rendering        - http://www.bonzaisoftware.com/tnp/gl-water-tutorial/#op
// Normal calc in shader                - http://stackoverflow.com/questions/5281261/generating-a-normal-map-from-a-height-map
// Normal calc heightfield              - http://archive.gamedev.net/archive/reference/programming/features/normalheightfield/index.html
// Normal calc 2:                       - http://www.flipcode.com/archives/Calculating_Vertex_Normals_for_Height_Maps.shtml
// Caustics from height fields:         - http://research.cs.tamu.edu/keyser/Papers/CausticsCGI09.pdf


// v0.0.0.8 - https://gist.github.com/roxlu/06e3c60a587344ad2079 - working diffuse lighting
// v0.0.0.7 - https://gist.github.com/roxlu/3439fd568e9aaeae8b79 - working animated normals
// v0.0.0.2 - https://gist.github.com/roxlu/f232fae0ad07ed8c8b89
 

#define ROXLU_USE_OPENGL
#define ROXLU_USE_MATH
#define ROXLU_USE_PNG
#include <tinylib.h>
#include <vector>
 
extern "C" {
#  include "mikktspace.h"
}

#define FIELD_N 128
#define FIELD_NN  (FIELD_N * FIELD_N)

struct WaterVertex {
  WaterVertex(){}
  vec3 pos;
  vec2 tex;
  vec3 norm;
};

struct WaterFace {
  WaterFace():a(0),b(0),c(0){}
  WaterFace(int a, int b, int c):a(a),b(b),c(c){ indices[0] = a; indices[1] = b; indices[2] = c;}
  int a,b,c;
  vec3 norm;
  int indices[3];
};

static const char* WATER_VS = ""
  "#version 150\n"
  "uniform float u_time;"
  "uniform mat4 u_pm;"
  "uniform mat4 u_vm;"
  "in vec4 a_pos;"
  "in vec2 a_tex;"
  "in vec3 a_norm;"
  "out vec2 v_tex;"
  "out vec3 v_norm;"
  "out vec3 v_view_pos;"
  "out vec3 v_world_pos;"
  "void main() {"
  "  gl_Position = u_pm * u_vm * a_pos;"
  "  v_tex = a_tex;"
  "  v_norm = a_norm;"
  "  v_world_pos = a_pos.xyz;"
  "  v_view_pos = vec3(u_vm * a_pos);"
  "}"
  "";

static const char* WATER_FS = ""
  "#version 150\n"
  "uniform mat4 u_vm;"
  "uniform sampler2D u_norm_tex;"
  "uniform sampler2D u_flow_tex;"
  "uniform sampler2D u_noise_tex;"
  "uniform sampler2D u_diffuse_tex;"
  "uniform samplerCube u_cube_tex;"
  "uniform float u_time;"
  "uniform float u_time0;"
  "uniform float u_time1;"
  "out vec4 fragcolor;"
  "in vec2 v_tex;"
  "in vec3 v_norm;"
  "in vec3 v_view_pos;"
  "in vec3 v_world_pos;"

  "void main() {"

  // flow map
  "  fragcolor = vec4(1.0, 0.0, 0.0, 1.0);"
  "  vec4 norm_color = texture(u_norm_tex, v_tex);"
  "  vec4 flow_color = texture(u_flow_tex, v_tex);"
  "  vec4 noise_color = texture(u_noise_tex, v_tex);"
  "  vec4 diffuse_color = texture(u_diffuse_tex, v_tex);"
  "  vec2 flow = flow_color.rg;"
  "  float cycle_offset = noise_color.r;"
  "  float phase0 = (cycle_offset * 0.1 + u_time0);"
  "  float phase1 = (cycle_offset * 0.1 + u_time1);"
  "  float tex_scale = 1.0;"
  "  vec2 texcoord0 = (v_tex * tex_scale) + (flow * phase0);"
  "  vec2 texcoord1 = (v_tex * tex_scale) + (flow * phase1);"
  "  vec3 normal0 = texture(u_norm_tex, texcoord0).rgb;"
  "  vec3 normal1 = texture(u_norm_tex, texcoord1).rgb;"
  "  float lerp = abs((0.5 - u_time0)) / 0.5;"
  "  vec3 offset = (mix(normal0, normal1, lerp));"


  "  offset = normalize(offset);" // do we need to normalize? Possibly because else it will be weighter differently when we add v_norm
  "  offset = normalize(v_norm + offset);"
  "  offset = -1.0 + 2.0 * offset;"


  "  vec3 diffuse0 = texture(u_diffuse_tex, texcoord0).rgb;"
  "  vec3 diffuse1 = texture(u_diffuse_tex, texcoord1).rgb;"
  "  diffuse_color.rgb = (mix(diffuse0, diffuse1, lerp));"

#if 0
  "  fragcolor.rgb = offset;"
  "  fragcolor.rgb = diffuse_color.rgb;"
  "  fragcolor.rgb = 0.5 + 0.5 * offset;"
  "  fragcolor.rgb = normalize(v_norm);"
#else

  // create the tangent space
  " vec3 tang = normalize(mat3(u_vm) * vec3(0.0, 0.0, 1.0));"
  " vec3 norm = normalize(mat3(u_vm) * offset);"
  " vec3 binorm = normalize(cross(tang, norm));"
  " mat3 to_tang = mat3("
  "     tang.x, binorm.x, norm.x, "
  "     tang.y, binorm.y, norm.y, "
  "     tang.z, binorm.z, norm.z);"
  " vec3 light = vec3(0.0, 100.0, 0.0);"
  " vec3 pos = mat3(u_vm) * v_world_pos;"
  " vec3 light_dir = normalize(to_tang * (light - pos));"
  " vec3 view_dir = to_tang * normalize(-pos);"
  " vec3 r = reflect(-light_dir, norm);"
  " float sdn = max(dot(light_dir, norm), 0.0);"
  " vec3 spec = vec3(0.0);"
  " if(sdn > 0.0) {"
  "    spec = vec3(1.0) * pow(max(dot(r, view_dir), 0.0), 0.3);"
  " }"
  " fragcolor.rgb = diffuse_color.rgb + 1.5 * sdn * vec3(1.0, 0.0, 0.0) + spec * 0.5;"
  //  " fragcolor.rgb = spec;"
  //  " fragcolor.rgb = vec3(1.0, 0.0, 0.0) * (sdn);"
  
  #if 0
  // convert everything to eye space
  "  vec3 norm_eye = mat3(u_vm) * offset;"
  "  vec3 pos_eye = mat3(u_vm) * v_world_pos;"
  "  vec3 light_world_a = vec3(10.0, -100.0, sin(u_time) * 10.0);"
  "  vec3 light_world_a_eye = mat3(u_vm) * light_world_a;"
  "  vec3 dir_to_light_a = normalize(light_world_a_eye - pos_eye);"
  "  float ndl_a = max(dot(dir_to_light_a, norm_eye), 0.0);"
  "  vec3 light_color = vec3(1.0, 0.0, 0.0) * ndl_a; "

  // specular
  "  vec3 dir_from_light_eye = normalize(pos_eye - light_world_a_eye);"
  "  vec3 reflected = normalize(reflect(dir_from_light_eye, norm_eye));"
  "  float sdn = max(dot(dir_from_light_eye, norm_eye), 0.0);"
  "  vec3 spec_color = vec3(0.0);"
  "  if(sdn > 0.0) { "
  "     spec_color = pow(max(dot(reflected, normalize(pos_eye)), 0.0), 4.0) * vec3(1.0);"
  "  }"

  // cubemap
  "  vec3 cam_pos_world = vec3(0.0, 10.0, 0.0);"
  "  vec3 dir_to_cam_world = normalize(cam_pos_world - v_world_pos);"
  "  vec3 refract_world = refract(-dir_to_cam_world, offset, 1.333);"
  "  vec3 refract_color = texture(u_cube_tex, refract_world).rgb;"

  "  vec3 shaded_color = diffuse_color.rgb + light_color * 0.5 + spec_color * 0.1;"

  "  fragcolor.rgb = shaded_color;"
  #endif
#if 0
  "  if(gl_FragCoord.y < 384) {"
  "    fragcolor.rgb = spec_color;"
  //  "    fragcolor.rgb = ndl_a * vec3(1.0);"
  // "    fragcolor.rgb = normalize(v_norm);"
  "  }"
#endif

  //"  fragcolor.r = v_world_pos.y;"
  // "  vec3 shaded_color = diffuse_color.rgb + light_color * 0.4 + spec_color * 0.6 + refract_color * 0.3;"
    "  fragcolor.rgb = vec3(1.0, 0.0, 0.0);"
  //" fragcolor.rgb = 0.5 + 0.5 * offset;"
  //"fragcolor.rgb = spec_color;"

#endif
  "}"
  "";


// -----------------------------------------------------------
int water_get_num_faces(const SMikkTSpaceContext* ctx);
inline int water_get_num_vertices_of_face(const SMikkTSpaceContext* ctx, int iface) { return 3; }
void water_get_position(const SMikkTSpaceContext* ctx, float pos[], const int iface, const int ivert);
void water_get_normal(const SMikkTSpaceContext* ctx, float norm[], const int iface, const int ivert);
void water_get_texcoord(const SMikkTSpaceContext* ctx, float tex[], const int iface, const int ivert);
void water_set_tspace(const SMikkTSpaceContext* ctx, const float tangent[], const float sign, const int iface, const int ivert);
// -----------------------------------------------------------


class Water {

 public:
  Water();
  bool setup(int w, int h);
  void update(float dt);
  void draw();

  void testFlowFieldForce();
  void addDrop(int col, int row);

 private:

  /* Water simulation */
  void createVertices();
  void updateVertices();
  void calculateNormals();
  GLuint createTexture(std::string filename);
  GLuint createCubeMap(std::string prefix);

  /* Height field */
  bool setupHeightField();
  void updateHeightField(float dt);
  void applyHeightFieldToWater();

  /* Tangent space */
  void calculateTangents();

public:
  
  /* Water simulation */
  std::vector<WaterVertex> tmp_vertices;
  std::vector<WaterVertex> vertices;
  std::vector<WaterFace> faces;
  std::vector<GLuint> indices;
  size_t bytes_allocated;

  mat4 pm;
  mat4 vm;
  GLuint vbo;
  //GLuint vbo_els;
  GLuint vao;
  GLuint vert;
  GLuint frag;
  GLuint prog;
  GLuint flow_tex;
  GLuint norm_tex;
  GLuint noise_tex;
  GLuint diffuse_tex;
  GLuint cube_tex;

  /* Height field */
  float u0[FIELD_NN];
  float u1[FIELD_NN];
  float v[FIELD_NN];

  /* Normal and tangent space. */
  SMikkTSpaceContext mikk_context;
  SMikkTSpaceInterface mikk_interface;
};

#endif