Separate classes into modules

10 files changed, 640 insertions, 0 deletions

diff --git a/boltzgas/__init__.py b/boltzgas/__init__.py
new file mode 100644
index 0000000..3b623dd
--- /dev/null
+++ b/boltzgas/__init__.py
@@ -0,0 +1,4 @@
+import boltzgas.initial_condition
+import boltzgas.visual
+
+from boltzgas.simulation import HardSphereSetup, HardSphereSimulation
diff --git a/boltzgas/initial_condition.py b/boltzgas/initial_condition.py
new file mode 100644
index 0000000..4d66263
--- /dev/null
+++ b/boltzgas/initial_condition.py
@@ -0,0 +1,15 @@
+import numpy as np
+
+def grid_of_random_velocity_particles(width, radius, u_scale):
+    np_position = np.ndarray((width**2, 2))
+    np_velocity = np.ndarray((width**2, 2))
+
+    grid = np.meshgrid(np.linspace(2*radius, 1-2*radius, width),
+                       np.linspace(2*radius, 1-2*radius, width))
+    np_position[:,0] = grid[0].flatten()
+    np_position[:,1] = grid[1].flatten()
+
+    np_velocity[:,0] = u_scale*(-0.5 + np.random.random_sample((width**2,)))
+    np_velocity[:,1] = u_scale*(-0.5 + np.random.random_sample((width**2,)))
+
+    return (np_position, np_velocity)
diff --git a/boltzgas/kernel.cl b/boltzgas/kernel.cl
new file mode 100644
index 0000000..ad984bb
--- /dev/null
+++ b/boltzgas/kernel.cl
@@ -0,0 +1,156 @@
+#define RADIUS      $radius
+#define N_PARTICLES $n_particles
+#define DELTA_T     $delta_t
+
+typedef float  scalar_t;
+typedef float2 vec_t;
+
+#define SCALAR_MAX FLT_MAX
+
+__constant vec_t wall_normals[4] = {
+	(vec_t)(-1.,0.),
+	(vec_t)( 1.,0.),
+	(vec_t)(0., 1.),
+	(vec_t)(0.,-1.),
+};
+
+__constant scalar_t wall_loc[4] = {
+	-RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS
+};
+
+scalar_t solve_wall_collision(vec_t n, scalar_t loc, vec_t p, vec_t v) {
+	if (dot(n,v) > 0.) {
+		vec_t wall_v = dot(n,v) * n;
+		return (loc - dot(p,n)) / dot(wall_v,n);
+	} else {
+		return SCALAR_MAX;
+	}
+}
+
+scalar_t pos_min_or_infty(scalar_t t0, scalar_t t1) {
+	if (t0 >= 0.) {
+		if (t1 >= 0.) {
+			return min(min(t0, t1), SCALAR_MAX);
+		} else {
+			return min(t0, SCALAR_MAX);
+		}
+	} else {
+		if (t1 >= 0.) {
+			return min(t1, SCALAR_MAX);
+		} else {
+			return SCALAR_MAX;
+		}
+	}
+}
+
+scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) {
+	scalar_t a = dot(v-v_, v-v_);
+	scalar_t b  = 2.*dot(p-p_, v-v_);
+	scalar_t c = dot(p-p_, p-p_) - 4.*RADIUS*RADIUS;
+	scalar_t d = b*b - 4.*a*c;
+
+	if (d >= 0.) {
+		scalar_t t0 = (-b + sqrt(d))/(2.*a);
+		scalar_t t1 = (-b - sqrt(d))/(2.*a);
+
+		return pos_min_or_infty(t0, t1);
+	} else {
+		return SCALAR_MAX;
+	}
+}
+
+__kernel void evolve(__global vec_t* pos_a,
+                     __global vec_t* vel_a,
+                     __global vec_t* pos_b,
+                     __global vec_t* vel_b,
+                     __global volatile unsigned int* last_collide)
+{
+	unsigned int i = get_global_id(0);
+
+	vec_t p = pos_a[i];
+	vec_t v = vel_a[i];
+
+	unsigned int jParticle = N_PARTICLES;
+	scalar_t min2intersect = SCALAR_MAX;
+
+	for (unsigned int iParticle=0; iParticle < N_PARTICLES; ++iParticle) {
+		if (iParticle != i && !(last_collide[i] == iParticle && last_collide[iParticle] == i)) {
+			vec_t p_ = pos_a[iParticle];
+			vec_t v_ = vel_a[iParticle];
+
+			scalar_t time2intersect = solve_collision(p, v, p_, v_);
+
+			if (time2intersect < min2intersect) {
+				min2intersect = time2intersect;
+				jParticle = iParticle;
+			}
+		}
+	}
+
+	unsigned int jWall = N_PARTICLES;
+	scalar_t min2wall = SCALAR_MAX;
+
+	for (unsigned int iWall=0; iWall < 4; ++iWall) {
+		scalar_t time2wall = solve_wall_collision(wall_normals[iWall], wall_loc[iWall], p, v);
+		if (time2wall < min2wall) {
+			min2wall = time2wall;
+			jWall = iWall;
+		}
+	}
+
+	if (min2intersect < DELTA_T) {
+		if (min2wall < min2intersect) {
+			p += min2wall * v;
+			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			p += (DELTA_T - min2wall) * v;
+			last_collide[i] = N_PARTICLES;
+
+			pos_b[i] = p;
+			vel_b[i] = v;
+		} else {
+			if (i < jParticle) {
+				vec_t p_ = pos_a[jParticle];
+				vec_t v_ = vel_a[jParticle];
+
+				p  += min2intersect * v;
+				p_ += min2intersect * v_;
+
+				vec_t omega = normalize(p - p_);
+
+				v  -= dot(vel_a[i] - vel_a[jParticle], omega) * omega;
+				v_ -= dot(vel_a[jParticle] - vel_a[i], omega) * omega;
+
+				p  += (DELTA_T - min2intersect) * v;
+				p_ += (DELTA_T - min2intersect) * v_;
+
+				pos_b[i] = p;
+				vel_b[i] = v;
+
+				pos_b[jParticle] = p_;
+				vel_b[jParticle] = v_;
+
+				last_collide[i] = jParticle;
+				last_collide[jParticle] = i;
+			}
+		}
+	} else {
+		if (min2wall < DELTA_T) {
+			p += min2wall * v;
+			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			p += (DELTA_T - min2wall) * v;
+			last_collide[i] = N_PARTICLES;
+		} else {
+			p += DELTA_T * v;
+		}
+
+		pos_b[i] = p;
+		vel_b[i] = v;
+	}
+}
+
+__kernel void get_velocity_norms(__global vec_t* velocities, __global scalar_t* norms)
+{
+	unsigned int i = get_global_id(0);
+
+	norms[i] = length(velocities[i]);
+}
diff --git a/boltzgas/simulation.py b/boltzgas/simulation.py
new file mode 100644
index 0000000..a4ae4da
--- /dev/null
+++ b/boltzgas/simulation.py
@@ -0,0 +1,131 @@
+import numpy as np
+
+import pyopencl as cl
+mf = cl.mem_flags
+from pyopencl.tools import get_gl_sharing_context_properties
+
+import OpenGL.GL as gl
+from OpenGL.arrays import vbo
+
+from string import Template
+
+class HardSphereSetup:
+    def __init__(self, radius, char_u, position, velocity):
+        self.radius = radius
+        self.char_u = char_u
+        self.position = position
+        self.velocity = velocity
+        self.n_particles = len(position[:,0])
+
+def build_kernel(delta_t, n_particles, radius):
+    with open('boltzgas/kernel.cl', 'r') as kernel_src:
+        return Template(kernel_src.read()).substitute(
+            delta_t     = delta_t,
+            n_particles = n_particles,
+            radius      = radius)
+
+class HardSphereSimulation:
+    def setup_cl(self):
+        self.platform = cl.get_platforms()[0]
+        if self.opengl:
+            self.context = cl.Context(
+                properties=[(cl.context_properties.PLATFORM, self.platform)] + get_gl_sharing_context_properties())
+        else:
+            self.context = cl.Context(
+                properties=[(cl.context_properties.PLATFORM, self.platform)])
+        self.queue = cl.CommandQueue(self.context)
+
+        self.program = cl.Program(self.context, self.kernel_src).build(
+            '-cl-single-precision-constant -cl-opt-disable')
+
+        self.cl_particle_velocity_a = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity)
+        self.cl_particle_velocity_b = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity)
+
+        if self.opengl:
+            self.gl_particle_position_a = vbo.VBO(data=self.np_particle_position, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
+            self.gl_particle_position_a.bind()
+            self.gl_particle_position_b = vbo.VBO(data=self.np_particle_position, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
+            self.gl_particle_position_b.bind()
+
+            self.cl_particle_position_a = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particle_position_a))
+            self.cl_particle_position_b = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particle_position_b))
+        else:
+            self.cl_particle_position_a = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_position)
+            self.cl_particle_position_b = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_position)
+
+        self.cl_last_collide = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_last_collide)
+        self.cl_particle_velocity_norms = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity_norms)
+
+    def __init__(self, setup, opengl = False, t_scale = 1.0):
+        self.np_particle_position = setup.position.astype(np.float32)
+        self.np_particle_velocity = setup.velocity.astype(np.float32)
+
+        self.n_particles = setup.n_particles
+        self.radius      = setup.radius
+        self.char_u      = setup.char_u
+
+        self.opengl = opengl
+        self.t_scale = t_scale
+
+        self.np_last_collide = np.ndarray((self.n_particles, 1), dtype=np.uint32)
+        self.np_last_collide[:,0] = self.n_particles
+
+        self.np_particle_velocity_norms = np.ndarray((self.n_particles, 1), dtype=np.float32)
+
+        self.kernel_src = build_kernel(self.t_scale*self.radius/self.char_u, self.n_particles, self.radius)
+
+        self.setup_cl()
+
+        self.tick = True
+
+    def evolve(self):
+        if self.opengl:
+            cl.enqueue_acquire_gl_objects(self.queue, [self.cl_particle_position_a, self.cl_particle_position_b])
+
+        if self.tick:
+            self.tick = False
+            kernelargs = (self.cl_particle_position_a, self.cl_particle_velocity_a, self.cl_particle_position_b, self.cl_particle_velocity_b, self.cl_last_collide)
+        else:
+            self.tick = True
+            kernelargs = (self.cl_particle_position_b, self.cl_particle_velocity_b, self.cl_particle_position_a, self.cl_particle_velocity_a, self.cl_last_collide)
+
+        self.program.evolve(self.queue, (self.n_particles,), None, *(kernelargs))
+
+    def gl_draw_particles(self):
+        gl.glEnableClientState(gl.GL_VERTEX_ARRAY)
+
+        if self.tick:
+            self.gl_particle_position_b.bind()
+            gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_b)
+        else:
+            self.gl_particle_position_a.bind()
+            gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_a)
+
+        gl.glDrawArrays(gl.GL_POINTS, 0, self.n_particles)
+        gl.glDisableClientState(gl.GL_VERTEX_ARRAY)
+
+    def get_velocities(self):
+        if self.tick:
+            cl.enqueue_copy(self.queue, self.np_particle_velocity, self.cl_particle_velocity_b).wait()
+            return self.np_particle_velocity
+        else:
+            cl.enqueue_copy(self.queue, self.np_particle_velocity, self.cl_particle_velocity_a).wait()
+            return self.np_particle_velocity
+
+    def get_velocity_norms(self):
+        if self.tick:
+            self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_b, self.cl_particle_velocity_norms)
+            cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait()
+            return self.np_particle_velocity_norms
+        else:
+            self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_a, self.cl_particle_velocity_norms)
+            cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait()
+            return self.np_particle_velocity_norms
+
+    def get_positions(self):
+        if self.tick:
+            cl.enqueue_copy(self.queue, self.np_particle_position, self.cl_particle_position_b).wait()
+            return self.np_particle_position
+        else:
+            cl.enqueue_copy(self.queue, self.np_particle_position, self.cl_particle_position_a).wait()
+            return self.np_particle_position
diff --git a/boltzgas/visual/__init__.py b/boltzgas/visual/__init__.py
new file mode 100644
index 0000000..3462d96
--- /dev/null
+++ b/boltzgas/visual/__init__.py
@@ -0,0 +1,4 @@
+from .box import ColoredBox
+from .histogram import VelocityHistogram
+from .tracer import Tracer
+from .view import View
diff --git a/boltzgas/visual/box.py b/boltzgas/visual/box.py
new file mode 100644
index 0000000..0a6ab5a
--- /dev/null
+++ b/boltzgas/visual/box.py
@@ -0,0 +1,16 @@
+from OpenGL.GL import *
+
+class ColoredBox:
+    def __init__(self, origin, extend, color):
+        self.origin = origin
+        self.extend = extend
+        self.color = color
+
+    def display(self, uniform):
+        glUniform3f(uniform['face_color'], *self.color)
+        glBegin(GL_TRIANGLE_STRIP)
+        glVertex(self.origin[0],                  self.origin[1]                 , 0.)
+        glVertex(self.origin[0] + self.extend[0], self.origin[1]                 , 0.)
+        glVertex(self.origin[0]                 , self.origin[1] + self.extend[1], 0.)
+        glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.)
+        glEnd()
diff --git a/boltzgas/visual/histogram.py b/boltzgas/visual/histogram.py
new file mode 100644
index 0000000..2f90b8b
--- /dev/null
+++ b/boltzgas/visual/histogram.py
@@ -0,0 +1,125 @@
+from OpenGL.GL import *
+
+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 concurrent.futures import Future, ProcessPoolExecutor
+
+def get_histogram(velocities, char_u):
+    maxwellian = stats.maxwell.fit(velocities)
+
+    plt.style.use('dark_background')
+
+    fig = plt.figure(figsize=(8,8))
+    ax = fig.add_axes([0.1, 0.06, 0.88, 0.92])
+
+    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()
+    width, height = fig.canvas.get_width_height()
+    texture = np.frombuffer(buf, dtype=np.uint8).reshape(width, height, 3)
+
+    plt.close(fig)
+
+    return (texture, height, width)
+
+
+class VelocityHistogram:
+    def __init__(self, gas, origin, extend):
+        self.gas = gas
+        self.origin = origin
+        self.extend = extend
+        self.steps = -1
+
+        self.pool = ProcessPoolExecutor(max_workers=1)
+        self.plotter = None
+
+        self.tick = False
+        self.mixing = 0.0
+
+    def setup(self):
+        self.vertices = np.array([
+            self.origin[0]                 , self.origin[1]                 , 0., 1.,
+            self.origin[0] + self.extend[0], self.origin[1]                 , 1., 1.,
+            self.origin[0]                 , self.origin[1] + self.extend[1], 0., 0.,
+            self.origin[0] + self.extend[0], self.origin[1] + self.extend[1], 1., 0.
+        ], dtype=np.float32)
+
+        self.vao = glGenVertexArrays(1)
+        self.vbo = glGenBuffers(1)
+
+        glBindVertexArray(self.vao)
+        glBindBuffer(GL_ARRAY_BUFFER, self.vbo)
+        glBufferData(GL_ARRAY_BUFFER, self.vertices.tostring(), GL_STATIC_DRAW)
+
+        glEnableVertexAttribArray(0)
+        glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(0))
+        glEnableVertexAttribArray(1)
+        glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(2*np.dtype(np.float32).itemsize))
+
+        self.texture_id = glGenTextures(2)
+
+        glBindTexture(GL_TEXTURE_2D, self.texture_id[0])
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+
+        glBindTexture(GL_TEXTURE_2D, self.texture_id[1])
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+
+    def update(self):
+        self.steps = self.steps + 1
+
+        if self.steps % 50 == 0 and self.plotter == None:
+            self.plotter = self.pool.submit(get_histogram, self.gas.get_velocity_norms(), self.gas.char_u)
+
+        else:
+            if self.plotter != None and self.plotter.done():
+                texture, width, height = self.plotter.result()
+                if self.tick:
+                    glBindTexture(GL_TEXTURE_2D, self.texture_id[0])
+                    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture)
+                    self.tick = False
+                    self.mixing = 1.0
+
+                else:
+                    glBindTexture(GL_TEXTURE_2D, self.texture_id[1])
+                    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture)
+                    self.tick = True
+                    self.mixing = 0.0
+
+                self.plotter = None
+
+
+    def display(self, uniform):
+        if self.tick:
+            self.mixing = min(self.mixing+0.1, 1.0);
+        else:
+            self.mixing = max(self.mixing-0.1, 0.0);
+
+        glBindTextures(self.texture_id[0], 2, self.texture_id)
+        glUniform1iv(uniform['picture'], len(self.texture_id), self.texture_id)
+        glUniform1f(uniform['mixing'], self.mixing)
+
+        glBindVertexArray(self.vao);
+        glDrawArrays(GL_TRIANGLE_STRIP, 0, 4)
+        glBindVertexArray(0)
+
diff --git a/boltzgas/visual/shader.py b/boltzgas/visual/shader.py
new file mode 100644
index 0000000..b4295a8
--- /dev/null
+++ b/boltzgas/visual/shader.py
@@ -0,0 +1,15 @@
+from OpenGL.GL import GL_VERTEX_SHADER, GL_FRAGMENT_SHADER
+from OpenGL.GL import shaders
+
+class Shader:
+    def __init__(self, fragment_src, vertex_src, uniform):
+        self.program = shaders.compileProgram(
+            shaders.compileShader(vertex_src, GL_VERTEX_SHADER),
+            shaders.compileShader(fragment_src, GL_FRAGMENT_SHADER))
+        self.uniform = { }
+        for name in uniform:
+            self.uniform[name] = shaders.glGetUniformLocation(self.program, name)
+
+    def use(self):
+        shaders.glUseProgram(self.program)
+
diff --git a/boltzgas/visual/tracer.py b/boltzgas/visual/tracer.py
new file mode 100644
index 0000000..a5b2203
--- /dev/null
+++ b/boltzgas/visual/tracer.py
@@ -0,0 +1,20 @@
+from OpenGL.GL import *
+
+class Tracer:
+    def __init__(self, gas, iParticle):
+        self.gas = gas
+        self.iParticle = iParticle
+        self.trace = [ ]
+
+    def update(self):
+        position = self.gas.get_positions()[self.iParticle]
+        self.trace.append((position[0], position[1]))
+
+    def display(self, uniform):
+        glUniform3f(uniform['face_color'], 1., 0., 0.)
+        glLineWidth(2)
+        glBegin(GL_LINE_STRIP)
+        for v in self.trace:
+            glVertex(*v, 0.)
+        glEnd()
+
diff --git a/boltzgas/visual/view.py b/boltzgas/visual/view.py
new file mode 100644
index 0000000..d6e188b
--- /dev/null
+++ b/boltzgas/visual/view.py
@@ -0,0 +1,154 @@
+from OpenGL.GL   import *
+from OpenGL.GLUT import *
+
+from pyrr import matrix44
+
+import numpy as np
+
+from .shader import Shader
+
+particle_shader = (
+    """
+        #version 430
+
+        in vec3 color;
+
+        void main(){
+            if (length(gl_PointCoord - vec2(0.5)) > 0.5) {
+                discard;
+            }
+
+            gl_FragColor = vec4(color.xyz, 0.0);
+        }
+    """,
+    """
+        #version 430
+
+        layout (location=0) in vec2 particles;
+
+        out vec3 color;
+
+        uniform mat4 projection;
+        uniform vec3 face_color;
+        uniform vec3 trace_color;
+        uniform uint trace_id;
+
+        void main() {
+            gl_Position = projection * vec4(particles, 0., 1.);
+
+            if (gl_VertexID == trace_id) {
+                color = trace_color;
+            } else {
+                color = face_color;
+            }
+        }
+    """,
+    ['projection', 'face_color', 'trace_color', 'trace_id']
+)
+
+decoration_shader = (
+    """
+        #version 430
+
+        in vec3 color;
+
+        void main(){
+            gl_FragColor = vec4(color.xyz, 0.0);
+        }
+    """,
+    """
+        #version 430
+
+        in vec3 vertex;
+
+        out vec3 color;
+
+        uniform mat4 projection;
+        uniform vec3 face_color;
+
+        void main() {
+            gl_Position = projection * vec4(vertex, 1.);
+            color = face_color;
+        }
+    """,
+    ['projection', 'face_color']
+)
+
+texture_shader = (
+    """
+        #version 430
+
+        in vec2 tex_coord;
+
+        uniform sampler2D picture[2];
+
+        uniform float mixing;
+
+        void main() {
+            gl_FragColor = mix(texture(picture[0], tex_coord), texture(picture[1], tex_coord), mixing);
+        }
+    """,
+    """
+        #version 430
+
+        layout (location=0) in vec2 screen_vertex;
+        layout (location=1) in vec2 texture_vertex;
+
+        out vec2 tex_coord;
+
+        uniform mat4 projection;
+
+        void main() {
+            gl_Position = projection * vec4(screen_vertex, 0.0, 1.0);
+            tex_coord = texture_vertex;
+        }
+    """,
+    ['picture','projection','mixing']
+)
+
+class View:
+    def __init__(self, gas, decorations, windows):
+        self.gas = gas
+        self.decorations = decorations
+        self.windows = windows
+
+        self.texture_shader    = Shader(*texture_shader)
+        self.particle_shader   = Shader(*particle_shader)
+        self.decoration_shader = Shader(*decoration_shader)
+
+    def reshape(self, width, height):
+        glViewport(0,0,width,height)
+
+        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.05, -1.0/2, 0])
+
+        self.pixels_per_unit = height / world_height
+        self.projection = np.matmul(translation, projection)
+
+    def display(self):
+        glClearColor(0.,0.,0.,1.)
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
+
+        self.decoration_shader.use()
+        glUniformMatrix4fv(self.decoration_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+        for decoration in self.decorations:
+            decoration.display(self.decoration_shader.uniform)
+
+        self.texture_shader.use()
+        glUniformMatrix4fv(self.texture_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+        for window in self.windows:
+            window.display(self.texture_shader.uniform)
+
+        self.particle_shader.use()
+        glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+        glUniform3f(self.particle_shader.uniform['face_color'],  1., 1., 1.)
+        glUniform3f(self.particle_shader.uniform['trace_color'], 1., 0., 0.)
+        glUniform1ui(self.particle_shader.uniform['trace_id'], 4)
+        glEnable(GL_POINT_SPRITE)
+        glPointSize(2*self.gas.radius*self.pixels_per_unit)
+        self.gas.gl_draw_particles()
+
+        glutSwapBuffers()