Plot and display histogram on texture

Hacked together but works well enough for demonstration

1 files changed, 126 insertions, 8 deletions

diff --git a/gas.py b/gas.py
index de047a1..c4c7062 100644
--- a/gas.py
+++ b/gas.py
@@ -6,8 +6,17 @@ from pyrr import matrix44
 
 import numpy as np
 
+import scipy.stats as stats
+import scipy.constants as const
+from scipy.optimize import minimize
+
+import matplotlib
+import matplotlib.pyplot as plt
+
 from particles import GasFlow, HardSphereSetup, grid_of_random_velocity_particles
 
+from concurrent.futures import Future, ThreadPoolExecutor
+
 glutInit()
 glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)
 glutInitWindowPosition(0, 0)
@@ -80,19 +89,43 @@ decoration_shader = Shader(
     uniform = ['projection', 'face_color']
 )
 
+texture_shader = Shader(
+    fragment_src = """
+        #version 130
+
+        uniform sampler2D screenTexture;
+
+        void main() {
+            gl_FragColor = texture(screenTexture, vec2(gl_TexCoord[0].x, 1.0 - gl_TexCoord[0].y));
+        }""",
+    vertex_src = """
+        #version 130
+
+        in vec3 vertex;
+
+        uniform mat4 projection;
+
+        void main() {
+            gl_Position = projection * vec4(vertex, 1.0);
+            gl_TexCoord[0] = gl_MultiTexCoord0;
+        }""",
+    uniform = ['projection']
+)
+
 class View:
-    def __init__(self, gas, decorations):
+    def __init__(self, gas, decorations, windows):
         self.gas = gas
         self.decorations = decorations
+        self.windows = windows
 
     def reshape(self, width, height):
         glViewport(0,0,width,height)
 
-        world_height = 1.0
+        world_height = 1.4
         world_width = world_height / height * width
 
         projection  = matrix44.create_orthogonal_projection(-world_width/2, world_width/2, -world_height/2, world_height/2, -1, 1)
-        translation = matrix44.create_from_translation([-1.0/2, -1.0/2, 0])
+        translation = matrix44.create_from_translation([-1.05, -1.0/2, 0])
 
         self.pixels_per_unit = height / world_height
         self.projection = np.matmul(translation, projection)
@@ -106,6 +139,11 @@ class View:
         for decoration in self.decorations:
             decoration.display(decoration_shader.uniform)
 
+        texture_shader.use()
+        glUniformMatrix4fv(texture_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+        for window in self.windows:
+            window.display(texture_shader.uniform)
+
         particle_shader.use()
         glUniformMatrix4fv(particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
         glEnable(GL_POINT_SPRITE)
@@ -147,26 +185,106 @@ class ColoredBox:
         glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.)
         glEnd()
 
+def get_histogram(velocities):
+    velocities = np.array([np.linalg.norm(v) for v in velocities])
+    maxwellian = stats.maxwell.fit(velocities)
+
+    fig = plt.figure(figsize=(10,10))
+
+    plt.ylim(0, 0.003)
+    plt.ylabel('Probability')
+
+    plt.xlim(0, 1.2*char_u)
+    plt.xlabel('Velocity magnitude [m/s]')
+
+    plt.hist(velocities, bins=50, density=True, alpha=0.5, label='Simulated velocities')
+
+    xs = np.linspace(0, 1.2*char_u, 100)
+    plt.plot(xs, stats.maxwell.pdf(xs, *maxwellian), label='Maxwell-Boltzmann distribution')
+
+    plt.legend(loc='upper right')
+
+    fig.canvas.draw()
+    buf = fig.canvas.tostring_rgb()
+    ncols, nrows = fig.canvas.get_width_height()
+    texture = np.frombuffer(buf, dtype=np.uint8).reshape(nrows, ncols, 3)
+    plt.close()
+
+    return (texture, ncols, nrows)
+
+class VelocityHistogram:
+    def __init__(self, gas, origin, extend):
+        self.gas = gas
+        self.origin = origin
+        self.extend = extend
+        self.steps = -1
+
+        self.pool = ThreadPoolExecutor(max_workers=1)
+        self.plotter = None
+
+    def setup(self):
+        self.texture_id = glGenTextures(1)
+        glBindTexture(GL_TEXTURE_2D, self.texture_id)
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+        glEnable(GL_TEXTURE_2D)
+
+    def update(self):
+        self.steps = self.steps + 1
+
+        if self.steps % 100 == 0 and self.plotter == None:
+            self.plotter = self.pool.submit(get_histogram, self.gas.get_velocities())
+
+        else:
+            if self.plotter != None:
+                if self.plotter.done():
+                    texture, width, height = self.plotter.result()
+                    self.plotter = None
+
+                    glEnable(GL_TEXTURE_2D)
+                    glBindTexture(GL_TEXTURE_2D, self.texture_id)
+                    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture)
+
+
+    def display(self, uniform):
+        glEnable(GL_TEXTURE_2D)
+        glBindTexture(GL_TEXTURE_2D, self.texture_id)
+
+        glBegin(GL_TRIANGLE_STRIP)
+        glTexCoord(0.,0.)
+        glVertex(self.origin[0],                  self.origin[1]                 , 0.)
+        glTexCoord(1.,0.)
+        glVertex(self.origin[0] + self.extend[0], self.origin[1]                 , 0.)
+        glTexCoord(0.,1.)
+        glVertex(self.origin[0]                 , self.origin[1] + self.extend[1], 0.)
+        glTexCoord(1.,1.)
+        glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.)
+        glEnd()
+
 grid_width = 30
-radius = 0.002
+radius = 0.003
 char_u = 1120
 
 position, velocity = grid_of_random_velocity_particles(grid_width, radius, char_u)
 velocity[:,:] = 0
 velocity[0,0] = 10*char_u
-velocity[0,1] =  1*char_u
+velocity[0,1] = 4*char_u
 
 config = HardSphereSetup(radius, char_u, position, velocity)
 gas = GasFlow(config, opengl = True, t_scale = 1.0)
 
 tracer = Tracer(gas, 5)
-view = View(gas, [ColoredBox([0,0], [1,1], (1,1,1)), tracer])
+histogram = VelocityHistogram(gas, [1.1,0], [1,1])
+histogram.setup()
+view = View(gas, [ ColoredBox([0,0], [1,1], (1,1,1)) ], [ histogram ])
 
 def on_display():
-    for i in range(0,10):
+    for i in range(0,2):
         gas.evolve()
 
-    tracer.update()
+    gas.queue.finish()
+
+    histogram.update()
 
     view.display()