diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0af8d15..74f67f4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.3.0] - 2022-09-23
+### Added
+- dsmc module
+
 ## [0.3.0] - 2022-09-22
 ### Added
 - octree module
diff --git a/dsmc/__init__.py b/dsmc/__init__.py
index e69de29..9ebec8a 100644
--- a/dsmc/__init__.py
+++ b/dsmc/__init__.py
@@ -0,0 +1 @@
+from .dsmc import DSMC
diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py
index b30dac7..5af06e9 100644
--- a/dsmc/dsmc.py
+++ b/dsmc/dsmc.py
@@ -1,7 +1,144 @@
 import math
-import numpy
+import numpy as np
 from numba import njit
+from . import particles as prt
+from . import octree as oc
 
 @njit
-def test():
-	return 1.0
+def _push(velocities, positions, dt):  
+    return positions + velocities*dt
+
+@njit 
+def _boundary(velocities, positions, domain):
+    for p in range(len(positions)):
+        while not oc._is_inside(positions[p], domain):
+            for i in range(3):
+                if positions[p][i] < domain[i][0]:
+                    positions[p][i] = 2.0 * domain[i][0] - positions[p][i]
+                    velocities[p][i] *= -1.0
+                if positions[p][i] > domain[i][1]:
+                    positions[p][i] = 2.0 * domain[i][1] - positions[p][i]
+                    velocities[p][i] *= -1.0
+                    
+    return (velocities, positions)
+ 
+@njit   
+def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single:
+    """
+    Parameters
+    ----------
+    vel1 : velocity
+    vel2 : velocity
+    sigma_T : float
+        total cross section
+    Vc : float
+        cell volume
+    w : float
+        weight
+    N : int
+        number of particles
+        
+    Returns
+    -------
+    collision proability : float
+    """
+    return np.linalg.norm(vel1 - vel2) * sigma_T * dt * w * N / Vc;
+  
+@njit  
+def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 : float, mass2 : float, rel_vel_module : float, rand_number1 : float, rand_number2 : float) -> tuple[np.ndarray, np.ndarray]:
+    mass12 = (mass1 + mass2)
+    mass1_12 = (mass1 / mass12)
+    mass2_12 = (mass2 / mass12)
+
+    cos_xi = (2.0 * rand_number1 - 1.0)
+    sin_xi = (np.sqrt(1.0 - cos_xi * cos_xi))
+    epsilon = (2.0 * np.pi * rand_number2)
+
+    centre_of_mass_velocity = (velocity1 * mass1 + velocity2 * mass2) * (1.0 / mass12)
+    
+    rel_velocity_new = np.empty((3, ))
+
+    rel_velocity_new[0] = rel_vel_module * cos_xi
+    rel_velocity_new[1] = rel_vel_module * sin_xi * np.cos(epsilon)
+    rel_velocity_new[2] = rel_vel_module * sin_xi * np.sin(epsilon)
+
+    return (centre_of_mass_velocity + rel_velocity_new * mass2_12 , centre_of_mass_velocity - rel_velocity_new * mass1_12)
+   
+@njit 
+def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, Vc : float, dt : float, w : float, offset : int, N : int) -> np.ndarray:
+    i = 1
+    while i < N:
+        p1 = permutations[offset + i - 1]
+        p2 = permutations[offset + i]
+        P = _calc_prob(velocities[p1], velocities[p2], sigma_T, Vc, dt, w, N)
+        R = np.random.random(3)
+        
+        if R[0] < P:
+            new_vels = _calc_post_col_vels(velocities[p1], velocities[p2], mass, mass, np.linalg.norm(velocities[p1] - velocities[p2]), R[1], R[2])
+            velocities[p1] = new_vels[0]
+            velocities[p2] = new_vels[1]
+        
+        i += 2
+    
+    return velocities
+ 
+@njit   
+def _update_vels(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, dt : float, w : float, elem_offsets : np.ndarray, number_elements : np.ndarray, number_children : np.ndarray, cell_boxes : np.ndarray, Nleafs : int) -> np.ndarray:
+    for i in range(Nleafs):
+        if not number_children[i]:
+            Vc = oc.get_V(cell_boxes[i])
+            velocities = _update_velocities(permutations, velocities, mass, sigma_T, Vc, dt, w, elem_offsets[i], number_elements[i])
+            
+    return velocities
+
+class DSMC:
+    def __init__(self):
+        self.clear()
+        
+    def clear(self):
+        self.particles = prt.Particles()
+        self.octree = oc.Octree()
+        self.w = None
+        self.domain = None
+        self.sigma_T = 3.631681e-19
+        self.mass = None
+        
+    def advance(self, dt):
+        if self.domain is None:
+            raise Exception("simulation domain not defined")
+        if self.particles.N == 0:
+            raise Exception("no particles created")
+        if self.w == None:
+            raise Exception("particle weight not set")
+            
+        self.octree.build(self.particles.Pos)
+        
+        self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos)
+        positions = _push(self.particles.Vel, self.particles.Pos, dt)
+        self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain)
+        
+    def _update_velocities(self, dt):
+        Nleafs : int = len(self.octree.leafs)
+        elem_offsets : np.ndarray = np.array([leaf.elem_offset for leaf in self.octree.leafs], dtype=int)
+        number_elements : np.ndarray = np.array([leaf.number_elements for leaf in self.octree.leafs], dtype=int)
+        number_children : np.ndarray = np.array([leaf.number_children for leaf in self.octree.leafs], dtype=int)
+        cell_boxes : np.ndarray = np.array([box for box in self.octree.cell_boxes])
+        
+        return _update_vels(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, dt, self.w, elem_offsets, number_elements, number_children, cell_boxes, Nleafs)
+        
+    def create_particles(self, box, T, n):
+        N = int(round(n / self.w))
+        print("creating {} particles".format(N))
+        self.particles.create_particles(box, self.mass, T, N)
+        
+        print("now containing {} particles, {} total".format(N, self.particles.N))
+        
+    def set_domain(self, domain):
+        self.domain = np.array(domain)
+        
+    def set_mass(self, mass):
+        self.mass = mass
+        
+    def set_weight(self, w):
+        self.octree.w = w
+        self.w = w
diff --git a/dsmc/octree.py b/dsmc/octree.py
index 81b18ac..09b3c33 100644
--- a/dsmc/octree.py
+++ b/dsmc/octree.py
@@ -110,6 +110,10 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray
             Nnew += 1
             
     return new_permutations, Nnew
+
+@njit
+def get_V(box):
+    return (box[0][1] - box[0][0]) * (box[1][1] - box[1][0]) * (box[2][1] - box[2][0])
     
 class Leaf:
     def __init__(self):
@@ -137,6 +141,7 @@ def clear(self):
         self.level = 0
         
     def build(self, positions):
+        self.clear()
         self._create_root(positions)
         self.permutations = np.array([i for i in range(len(positions))])
         
diff --git a/dsmc/particles.py b/dsmc/particles.py
index 852ff0c..9ca72eb 100644
--- a/dsmc/particles.py
+++ b/dsmc/particles.py
@@ -121,14 +121,18 @@ def VelPos(self):
 
     @VelPos.setter
     def VelPos(self, vel_pos):
-        assert vel_pos[0] == vel_pos[1]
-        assert isinstance(vel_pos[0], np.array)
-        assert isinstance(vel_pos[1], np.array)
+        assert len(vel_pos[0]) == len(vel_pos[1])
+
         self._velocities = vel_pos[0]
         self._positions = vel_pos[1]
         self._N = len(self._positions)
-
-    def create_particles(self, X, mass, T, N):
-        self._velocities = get_velocities(T, mass, N)
-        self._positions = calc_positions(X[0], X[1], X[2], N)
-        self._N  = N
+        
+    def create_particles(self, X, mass, T, N):    
+        if self._N == 0:
+            self._velocities = get_velocities(T, mass, N)
+            self._positions = calc_positions(X[0], X[1], X[2], N)
+            self._N  = N
+        else:
+            self._velocities = np.concatenate((self._velocities, get_velocities(T, mass, N)))
+            self._positions = np.concatenate((self._positions, calc_positions(X[0], X[1], X[2], N)))
+            self._N  += N
diff --git a/examples/plot.py b/examples/plot.py
new file mode 100644
index 0000000..b4e91ae
--- /dev/null
+++ b/examples/plot.py
@@ -0,0 +1,38 @@
+import matplotlib.pyplot as plt
+import csv
+import numpy as np
+
+if __name__ == '__main__':
+    res = []
+
+    with open('test.csv') as file:
+        reader = csv.reader(file, delimiter=',')
+        res = []
+        
+        for line in reader:
+            l = [m for m in line if m]
+            data = (np.array(l)).astype(float)
+            data = np.sort(data)
+            
+            N = 100
+            sor = np.zeros((N, ))
+            x = np.zeros((N, ))
+            dx = 0.05/N
+            x[0] = dx
+            q = 0
+            
+            for i in range(len(data)):
+                while data[i] > x[q]:
+                    q += 1
+                    x[q] = x[q - 1] + dx
+                    
+                sor[q] += 1
+                    
+            
+            x.resize(q)
+            sor.resize(q)
+            
+            res.append((x, sor))
+            
+    plt.plot(res[-1][0], res[-1][1])
+    plt.show()
diff --git a/examples/shock_tube.py b/examples/shock_tube.py
new file mode 100644
index 0000000..9b61c08
--- /dev/null
+++ b/examples/shock_tube.py
@@ -0,0 +1,40 @@
+import dsmc
+
+if __name__ == '__main__':
+    # general parameters
+    solver = dsmc.DSMC()
+    domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3))
+    dt = 1e-7
+    w = 2.4134e+15
+    mass = 6.6422e-26
+    niter = 300
+    
+    # low denisty particles
+    nhigh = 2.5e+20
+    Thigh = 300
+    Boxhigh = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3))
+    
+    # high denisty particles
+    nlow = 2.5e+19
+    Tlow = 300
+    Boxlow = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 25e-3))
+    
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+    
+    solver.create_particles(Boxlow, Tlow, nlow)
+    solver.create_particles(Boxhigh, Thigh, nhigh)
+    
+    with open("test.csv", "w") as file:
+        for it in range(niter):
+            print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
+            solver.advance(dt)
+            
+            for pos in solver.particles.Pos:
+                file.write("{:.4e},".format(pos[2]))
+                
+            file.write("\n")
+
+    print("")
+    print('done')
diff --git a/setup.py b/setup.py
index 13e9750..7963a05 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import setup
 setup(
     name='dsmc',
-    version='0.3.0',
+    version='0.4.0',
     author='Leo Basov',
     python_requires='>=3.6, <4',
     packages=["dsmc"],
diff --git a/tests/unit/test_dsmc/dsmc.py b/tests/unit/test_dsmc/dsmc.py
index 128cccb..90cdf10 100644
--- a/tests/unit/test_dsmc/dsmc.py
+++ b/tests/unit/test_dsmc/dsmc.py
@@ -1,6 +1,39 @@
 import unittest
-from dsmc import dsmc
+import numpy as np
+from dsmc import dsmc as ds
 
 class TestDSMC(unittest.TestCase):
-    def test_test(self):
-        self.assertEqual(1.0, dsmc.test())
+    def test_Constructor(self):
+        ds.DSMC()
+        
+    def test__calc_prob(self):
+        vel1 : np.ndarray = np.array([1.0, 2.0, 3.0])
+        vel2 : np.ndarray = np.array([1.0, 2.0, 3.0])
+        sigma_T : float = 1.0
+        Vc : float = 1.0
+        dt : float = 1.0
+        w : float = 1.0
+        N : int = 1
+        
+        res = ds._calc_prob(vel1, vel2, sigma_T, Vc, dt, w, N)
+        
+        self.assertEqual(np.linalg.norm(vel1 - vel2), res)
+        
+    def test__calc_post_col_vels(self):
+        velocity1 : np.ndarray = np.array([1.0, 2.0, 3.0])
+        velocity2 : np.ndarray = np.array([1.0, 2.0, 3.0])
+        mass1 : float = 1.0
+        mass2 : float = 1.0
+        rel_vel_module : float = 1.0
+        rand_number1 : float = 1.0
+        rand_number2 : float = 1.0
+        
+        res = ds._calc_post_col_vels(velocity1, velocity2, mass1, mass2, rel_vel_module, rand_number1, rand_number2)
+        
+        self.assertEqual(1.5, res[0][0])
+        self.assertEqual(2.0, res[0][1])
+        self.assertEqual(3.0, res[0][2])
+        
+        self.assertEqual(0.5, res[1][0])
+        self.assertEqual(2.0, res[1][1])
+        self.assertEqual(3.0, res[1][2])
diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py
index e27497f..1921691 100644
--- a/tests/unit/test_dsmc/particles.py
+++ b/tests/unit/test_dsmc/particles.py
@@ -63,3 +63,14 @@ def test_create_particles(self):
             self.assertTrue(particles.Pos[i][0] >= x[0] and particles.Pos[i][0] <= x[1])
             self.assertTrue(particles.Pos[i][1] >= y[0] and particles.Pos[i][1] <= y[1])
             self.assertTrue(particles.Pos[i][2] >= z[0] and particles.Pos[i][2] <= z[1])
+            
+        particles.create_particles(X, mass, T, N)
+
+        self.assertEqual(2*N, len(particles.Pos))
+        self.assertEqual(2*N, len(particles.Vel))
+        self.assertEqual(2*N, particles.N)
+        
+        for i in range(2*N):
+            self.assertTrue(particles.Pos[i][0] >= x[0] and particles.Pos[i][0] <= x[1])
+            self.assertTrue(particles.Pos[i][1] >= y[0] and particles.Pos[i][1] <= y[1])
+            self.assertTrue(particles.Pos[i][2] >= z[0] and particles.Pos[i][2] <= z[1])