3 files changed, 47 insertions, 10 deletions

diff --git a/channel_3d_volumetric_rendering_gl_interop.py b/channel_3d_volumetric_rendering_gl_interop.py
index efef0a2..91a470e 100644
--- a/channel_3d_volumetric_rendering_gl_interop.py
+++ b/channel_3d_volumetric_rendering_gl_interop.py
@@ -150,24 +150,63 @@ vec3 palette(float x) {
     );
 }
 
+vec3 blueBlackRedPalette(float x) {
+	if ( x < 0.5 ) {
+		return mix(
+			vec3(0.0, 0.0, 1.0),
+			vec3(0.0, 0.0, 0.0),
+			2*x
+		);
+	} else {
+		return mix(
+			vec3(0.0, 0.0, 0.0),
+			vec3(1.0, 0.0, 0.0),
+			2*(x - 0.5)
+		);
+	}
+}
+
+vec3 v(float x, float y, float z) {
+    return texture(moments, unit(vec3(x,y,z))).yzw;
+}
+
+vec3 curl(vec3 p) {
+    const float h = 1./$size_x;
+
+    const float dyvz = (v(p.x,p.y+1,p.z).z - v(p.x,p.y-1,p.z).z) / (2.0*h);
+    const float dzvy = (v(p.x,p.y,p.z+1).y - v(p.x,p.y,p.z-1).y) / (2.0*h);
+
+    const float dzvx = (v(p.x,p.y,p.z+1).x - v(p.x,p.y,p.z-1).x) / (2.0*h);
+    const float dxvz = (v(p.x+1,p.y,p.z).z - v(p.x-1,p.y,p.z).z) / (2.0*h);
+
+    const float dxvy = (v(p.x+1,p.y,p.z).y - v(p.x-1,p.y,p.z).y) / (2.0*h);
+    const float dyvx = (v(p.x,p.y+1,p.z).x - v(p.x,p.y-1,p.z).x) / (2.0*h);
+
+    return vec3(
+        dyvz - dzvy,
+        dzvx - dxvz,
+        dxvy - dyvx
+    );
+}
+
 vec3 trace(vec3 pos, vec3 ray) {
     const float ray_length = $max_ray_length;
     const float delta      = 1./ray_length;
 
     vec3 color  = vec3(0.0);
-    float value = 0.0;
 
     for (float t = 0.0; t < 1.0; t += delta) {
         const vec3 sample_pos = unit(pos + t*ray_length*ray);
         if (length(sample_pos) < sqrt(3.1)) {
             const vec4 data = texture(moments, sample_pos);
             if (data[3] != 1.0) {
-                const float norm = length(data.yzw) / $inflow;
-                value += delta * 0.5 * norm;
+                const vec3 c = curl(pos + t*ray_length*ray);
+
+                color += delta * blueBlackRedPalette(0.5+clamp(c.y,-1.0,1.0));
             } else {
-                const vec3 n = normalize(-0.5 + data.xyz); // recover surface normal
+                const vec3 n = data.xyz; // recover surface normal
                 const float brightness = clamp(dot(n, ray), 0, 1);
-                color = vec3(max(0.3,brightness));
+                color += vec3(max(0.3,brightness));
                 break;
             }
         } else {
@@ -175,8 +214,6 @@ vec3 trace(vec3 pos, vec3 ray) {
         }
     }
 
-    color += palette(value);
-
     return color;
 }
 
diff --git a/template/opengl.mako b/template/opengl.mako
index 50cdbbc..52699d5 100644
--- a/template/opengl.mako
+++ b/template/opengl.mako
@@ -115,7 +115,7 @@ __kernel void collect_gl_moments_to_texture(__global __read_only  ${float_type}*
       if (material_up    != 5) { n.z = -1; }
       if (material_down  != 5) { n.z =  1; }
 
-      data.xyz = 0.5 + 0.5*n; // pack surface normal into texture
+      data.xyz = normalize(n);
       data.w   = 1.0; // signal impermeable material to raytracer
     }
 
diff --git a/utility/opengl.py b/utility/opengl.py
index afbfa00..ff88b17 100644
--- a/utility/opengl.py
+++ b/utility/opengl.py
@@ -75,7 +75,7 @@ class MomentsTexture:
         glBindTexture(self.gl_texture_type, self.gl_moments)
 
         if self.gl_texture_type == GL_TEXTURE_3D:
-            glTexImage3D(self.gl_texture_type, 0, GL_RGBA, self.lattice.memory.size_x, self.lattice.memory.size_y, self.lattice.memory.size_z, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
+            glTexImage3D(self.gl_texture_type, 0, GL_RGBA32F, self.lattice.memory.size_x, self.lattice.memory.size_y, self.lattice.memory.size_z, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE)
@@ -83,7 +83,7 @@ class MomentsTexture:
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_R,     GL_CLAMP_TO_EDGE)
             self.cl_gl_moments  = cl.GLTexture(self.lattice.context, mf.READ_WRITE, self.gl_texture_type, 0, self.gl_moments, 3)
         elif self.gl_texture_type == GL_TEXTURE_2D:
-            glTexImage2D(self.gl_texture_type, 0, GL_RGBA, self.lattice.memory.size_x, self.lattice.memory.size_y, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
+            glTexImage2D(self.gl_texture_type, 0, GL_RGBA32F, self.lattice.memory.size_x, self.lattice.memory.size_y, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
             glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE)