diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4772e1f..0af8d15 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-22
+### Added
+- octree module
+
 ## [0.2.0] - 2022-09-15
 ### Added
 - particle module
diff --git a/doc/octree_resolution.ipynb b/doc/octree_resolution.ipynb
new file mode 100644
index 0000000..a4cfe8b
--- /dev/null
+++ b/doc/octree_resolution.ipynb
@@ -0,0 +1,68 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "9db716f6-c301-4d2e-b06f-19cd98de2e5f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "772d42fa-dc37-48e9-a599-392077582e98",
+   "metadata": {},
+   "source": [
+    "The mean free path can be written as\n",
+    "$$\n",
+    "\\lambda = \\frac{1}{\\sqrt{2} \\cdot n_{c} \\cdot \\sigma_T}\n",
+    "$$\n",
+    "where $n_c$ is the number density in the cell and $\\sigma_T$ is the total cross section.\n",
+    "The number of particles in cell $N_{c}$ can be calculated as (asuming cubic cells)\n",
+    "$$\n",
+    "N_{c} = \\frac{1}{w} \\cdot L_c^3 \\cdot n_c\n",
+    "$$\n",
+    "where $L_C$ is the side length of the cell and $w$ being the particle weight.\n",
+    "Assuming that the cell side fullfills the criteria\n",
+    "$$\n",
+    "\\frac{L_c}{\\lambda} \\leq \\frac{1}{2}.\n",
+    "$$\n",
+    "This leads to the criterium\n",
+    "$$\n",
+    "N \\leq \\frac{\\sqrt{2}}{32 \\cdot w \\cdot \\sigma^2_T \\cdot n^3_c}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f03b0f7-64a9-48e5-989c-2e253ce09c1d",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/dsmc/octree.py b/dsmc/octree.py
new file mode 100644
index 0000000..81b18ac
--- /dev/null
+++ b/dsmc/octree.py
@@ -0,0 +1,207 @@
+import math
+import numpy as np
+import numpy.typing as npt
+from numba import njit
+
+fmin = np.finfo(float).min
+fmax = np.finfo(float).max
+
+@njit
+def _find_bounding_box(positions : npt.NDArray) -> npt.NDArray:
+    box = np.array([[fmax, fmin], [fmax, fmin], [fmax, fmin]])
+    
+    for pos in positions:
+        for i in range(3):
+            if pos[i] < box[i][0]:
+                box[i][0] = pos[i]
+            if pos[i] > box[i][1]:
+                box[i][1] = pos[i]
+    
+    return box
+
+@njit    
+def _calc_N_res(w : float, sigma_T : float, n : float) -> int:
+    """
+    Parameters
+    ----------
+    w : float
+        particle weight
+    sigma_t : float
+              total cross section [m^2]
+    n : float
+        number density [1/m^3]
+    """
+
+    return int(round(np.sqrt(2.0) / (32.0 * w * sigma_T**3 * n**2)))
+
+@njit  
+def _calc_n(box : npt.NDArray, N : float, w : float) -> float:
+    """Calculates number density in cell
+    
+    Parameters
+    ----------
+    box : np.array(3, 3)
+          cell
+    N : int
+        number of particles in cell
+    w : float
+        particle weight
+        
+    Returns
+    -------
+    number density : float
+    """
+    return np.prod(np.array([box[i][1] - box[i][0] for i in range(3)])) * N / w
+
+@njit  
+def _is_resolved(box : npt.NDArray, N : int, w : float, sigma_T : float, Nmin : int, Nmax : int) -> bool:
+    if N == 0:
+        return False
+
+    n = _calc_n(box, N, w)
+    Nres = _calc_N_res(w, sigma_T, n)
+    
+    return N > 2 * min(Nmin, max(Nmin, Nres))
+
+@njit
+def _is_inside(position : npt.NDArray, box : npt.NDArray) -> bool:
+    a : bool = position[0] > box[0][0] and position[0] <= box[0][1]
+    b : bool = position[1] > box[1][0] and position[1] <= box[1][1]
+    c : bool = position[2] > box[2][0] and position[2] <= box[2][1]
+
+    return a and b and c
+
+@njit    
+def _swap(arr, pos1, pos2):
+    temp = arr[pos2]
+    arr[pos2] = arr[pos1]
+    arr[pos1] = temp
+    
+@njit
+def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray, offset : int, N : int) -> tuple[npt.NDArray, int]:
+    '''sort particles in cell
+    
+    Parameters
+    ----------
+    box : np.ndarray
+          cell
+    positions : ndarray((3, ))
+        particle positions
+    offset : int
+        number offset
+    N : int
+        number of particles to be considered
+        
+    Returns
+    -------
+    new_permutations : ndarray[int]
+    N : int
+        number of found positions
+    '''
+    new_permutations = np.copy(permutations)
+    temp = np.empty((3,))
+    runner = offset
+    Nnew = 0
+    for i in range(offset, offset + N):
+        p = permutations[i]
+        if _is_inside(positions[p], box):
+            _swap(new_permutations, i, runner)
+            runner += 1
+            Nnew += 1
+            
+    return new_permutations, Nnew
+    
+class Leaf:
+    def __init__(self):
+        self.level = 0
+        self.elem_offset = 0
+        self.number_elements = 0
+        self.id_parent = None
+        self.id_first_child = None
+        self.number_children = 0
+    
+class Octree:
+    def __init__(self):
+        self.clear()
+        
+    def clear(self):
+        self.cell_boxes = []
+        self.leafs = []
+        self.sigma_T = 3.631681e-19
+        self.w = 1.0
+        self.Nmin = 8
+        self.Nmax = 64
+        self.max_level = 10
+        self.permutations = []
+        self.cell_offsets = []
+        self.level = 0
+        
+    def build(self, positions):
+        self._create_root(positions)
+        self.permutations = np.array([i for i in range(len(positions))])
+        
+        for level in range(self.max_level):
+            self.level += 1
+            self.cell_offsets.append(self.cell_offsets[-1])
+            for i in range(self.cell_offsets[level], self.cell_offsets[level + 1]):
+                self._progress(i, positions)
+   
+            if self.cell_offsets[level + 1] == self.cell_offsets[level + 2]:
+                break            
+        
+    def _create_root(self, positions):
+        box = _find_bounding_box(positions)
+        leaf = Leaf()
+        leaf.number_elements = len(positions)
+        
+        self.cell_offsets += [0, 1]
+        self.leafs.append(leaf)
+        self.cell_boxes.append(box)
+        
+    def _progress(self, leaf_id, positions):
+        leaf = self.leafs[leaf_id]
+        if _is_resolved(self.cell_boxes[leaf_id], leaf.number_elements, self.w, self.sigma_T, self.Nmin, self.Nmax):
+            leaf.number_children = 8
+            leaf.id_first_child = leaf_id + 1
+            self.cell_offsets[-1] += 8
+            self._add_boxes(self.cell_boxes[leaf_id])
+        else:
+            pass
+            
+        offset = 0
+            
+        for i in range(leaf.number_children):
+            new_leaf = Leaf()
+            new_leaf.level = leaf.level + 1
+            new_leaf.id_parent = leaf_id
+
+            self.permutations, N = _sort(self.permutations, self.cell_boxes[leaf_id + 1 + i], positions, leaf.elem_offset, leaf.number_elements)
+            
+            new_leaf.number_elements = N
+            new_leaf.elem_offset = leaf.elem_offset + offset
+            offset += N
+
+            self.leafs.append(new_leaf)
+       
+    def _add_boxes(self, box):
+        half = np.array([0.5*(box[i][0] + box[i][1]) for i in range(3)])
+        
+        child_geo1 = np.array(((half[0], box[0][1]), (half[1], box[1][1]), (half[2], box[2][1])))
+        child_geo2 = np.array(((box[0][0], half[0]), (half[1], box[1][1]), (half[2], box[2][1])))
+        child_geo3 = np.array(((box[0][0], half[0]), (box[1][0], half[1]), (half[2], box[2][1])))
+        child_geo4 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (half[2], box[2][1])))
+        
+        child_geo5 = np.array(((half[0], box[0][1]), (half[1], box[1][1]), (box[2][0], half[2])))
+        child_geo6 = np.array(((box[0][0], half[0]), (half[1], box[1][1]), (box[2][0], half[2])))
+        child_geo7 = np.array(((box[0][0], half[0]), (box[1][0], half[1]), (box[2][0], half[2])))
+        child_geo8 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (box[2][0], half[2])))
+        
+        self.cell_boxes.append(child_geo1)
+        self.cell_boxes.append(child_geo2)
+        self.cell_boxes.append(child_geo3)
+        self.cell_boxes.append(child_geo4)
+        
+        self.cell_boxes.append(child_geo5)
+        self.cell_boxes.append(child_geo6)
+        self.cell_boxes.append(child_geo7)
+        self.cell_boxes.append(child_geo8)
diff --git a/setup.py b/setup.py
index 6d1b41d..13e9750 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import setup
 setup(
     name='dsmc',
-    version='0.1.0',
+    version='0.3.0',
     author='Leo Basov',
     python_requires='>=3.6, <4',
     packages=["dsmc"],
diff --git a/tests/unit/main.py b/tests/unit/main.py
index d5c2891..da2edd5 100644
--- a/tests/unit/main.py
+++ b/tests/unit/main.py
@@ -4,6 +4,7 @@
 
 from test_dsmc.dsmc import *
 from test_dsmc.particles import *
+from test_dsmc.octree import *
 
 if __name__ == '__main__':
 	unittest.main()
diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py
new file mode 100644
index 0000000..1ad85f1
--- /dev/null
+++ b/tests/unit/test_dsmc/octree.py
@@ -0,0 +1,104 @@
+import unittest
+import numpy as np
+from dsmc import octree as oc
+from dsmc import particles as part
+
+class TestOctree(unittest.TestCase):
+    def test__find_bounding_box(self):
+        positions = np.array([(0.0, 0.0, -1.0), (-2.0, -3.0, 0.0), (4.0, 5.0, 0.0)])
+        box = oc._find_bounding_box(positions)
+        
+        self.assertEqual(-2.0, box[0][0])
+        self.assertEqual(4.0, box[0][1])
+        
+        self.assertEqual(-3.0, box[1][0])
+        self.assertEqual(5.0, box[1][1])
+        
+        self.assertEqual(-1.0, box[2][0])
+        self.assertEqual(0.0, box[2][1])
+        
+    def test__calc_N_res(self):
+        w = 1.0e+9
+        sigma_T = 1.0e-16
+        n = 1.0e+17
+        ref = int(round(np.sqrt(2) / (32.0 * w * sigma_T**3 * n**2)))
+        N = oc._calc_N_res(w, sigma_T, n)
+        
+        self.assertEqual(ref, N)
+        
+    def test__calc_n(self):
+        box = ((0, 1), (2, 4), (4, 7))
+        N = 300
+        w = 100
+        res = oc._calc_n(box, N, w)
+        
+        self.assertEqual(18.0, res)
+        
+    def test__is_inside(self):
+        box = [[2.0, 3.0], [4.0, 6.0], [-1.0, 1.0]]
+        box = np.array(box)
+        position1 = np.array([2.5, 5.0, 0.0])
+        position2 = np.array([0.0, 0.0, 0.0])
+        position3 = np.array([2.5, 6.0, 0.5])
+        
+        self.assertTrue(oc._is_inside(position1, box))
+        self.assertFalse(oc._is_inside(position2, box))
+        self.assertTrue(oc._is_inside(position3, box))
+        
+    def test_sort(self):
+        box = np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]])
+        positions1 = np.random.random((100, 3))
+        positions2 = np.random.random((200, 3)) - np.ones((200, 3))*2
+        positions = np.concatenate((positions2, positions1, positions2))
+        permutations1 = np.array([i for i in range(len(positions2))])
+        permutations2 = np.array([i for i in range(len(positions2), len(positions))])
+        np.random.shuffle(permutations2)
+        permutations = np.concatenate((permutations1, permutations2))
+        offset = len(positions2)
+        N = len(positions1) + len(positions2)
+        count = 0
+        
+        for position in positions1:
+            self.assertTrue(oc._is_inside(position, box))
+            
+        for position in positions2:
+            self.assertFalse(oc._is_inside(position, box))
+            
+        for position in positions:
+            if oc._is_inside(position, box):
+                count += 1
+                
+        self.assertEqual(len(positions1), count)
+        
+        self.assertEqual(len(positions2), len(permutations1))
+        self.assertEqual(len(positions1) + len(positions2), len(permutations2))
+        self.assertEqual(len(positions), len(permutations))
+        
+        permutations, Nnew = oc._sort(permutations, box, positions, offset, N)
+        
+        self.assertEqual(Nnew, len(positions1))
+        
+        for i in range(offset, offset + Nnew):
+            p = permutations[i]
+            self.assertTrue(oc._is_inside(positions[p], box))
+            
+        for i in range(offset + Nnew, len(permutations)):
+            p = permutations[i]
+            self.assertFalse(oc._is_inside(positions[p], box))
+            
+            
+class TestOctreeOctree(unittest.TestCase):
+    def test_build(self):
+        positions = np.random.random((1000, 3))*2.0 - np.ones((1000, 3))
+        octree = oc.Octree()
+        octree.w = 1e+18
+        
+        octree.build(positions)
+        
+    def test__add_boxes(self):
+        octree = oc.Octree()
+        box = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]])
+        
+        octree._add_boxes(box)
+        
+        self.assertEqual(8, len(octree.cell_boxes))