diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6236035..c85a38a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,9 +6,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-## [0.5.1] - 2022-09-24
+## [0.6.0] - 2022-09-28
 ### Added
-- fixed bugs in octree module
+- functionality for aspect ratio check in octree
+
+### Fixed
+- bug in temperature initialization
+- several bugs in octree
+
+## [0.5.1] - 2022-09-24
+### Fixed
+- bugs in octree module
 
 ## [0.5.0] - 2022-09-23
 ### Added
diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py
index dd6aeb4..56b0b7f 100644
--- a/dsmc/diagnostics.py
+++ b/dsmc/diagnostics.py
@@ -1,5 +1,4 @@
 import numpy as np
-from numba import njit
 from . import particles as prt
 from . import octree as oc
 
@@ -26,7 +25,6 @@ def sort_bin(positions, axis, Nbin):
     
 def calc_n(bins, box, axis, w):
     Nbins = len(bins)
-    dx = (box[axis][1] - box[axis][0]) / Nbins
     V = 1
     n = np.empty((Nbins, ))
     for i in range(3):
diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py
index ad25fa4..3c13811 100644
--- a/dsmc/dsmc.py
+++ b/dsmc/dsmc.py
@@ -1,16 +1,18 @@
-import math
 import numpy as np
 from numba import njit
+from numba import prange
 from . import particles as prt
 from . import octree as oc
 
-@njit
-def _push(velocities, positions, dt):  
-    return positions + velocities*dt
+@njit(parallel=True)
+def _push(velocities, positions, dt):
+    for p in prange(len(positions)):
+        positions[p] = positions[p] + velocities[p]*dt
+    return positions
 
-@njit 
+@njit(parallel=True)
 def _boundary(velocities, positions, domain):
-    for p in range(len(positions)):
+    for p in prange(len(positions)):
         while not oc._is_inside(positions[p], domain):
             for i in range(3):
                 if positions[p][i] < domain[i][0]:
@@ -19,11 +21,11 @@ def _boundary(velocities, positions, domain):
                 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:
+
+@njit
+def _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single:
     """
     Parameters
     ----------
@@ -37,14 +39,14 @@ def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float
         weight
     N : int
         number of particles
-        
+
     Returns
     -------
     collision proability : float
     """
-    return np.linalg.norm(vel1 - vel2) * sigma_T * dt * w * N / Vc;
-  
-@njit  
+    return rel_vel * 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)
@@ -55,7 +57,7 @@ def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 :
     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
@@ -63,38 +65,36 @@ def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 :
     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 
+
+@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:
+    for i in range(1, N, 2):
         p1 = permutations[offset + i - 1]
         p2 = permutations[offset + i]
-        P = _calc_prob(velocities[p1], velocities[p2], sigma_T, Vc, dt, w, N)
+        rel_vel = np.linalg.norm(velocities[p1] - velocities[p2])
+        P = _calc_prob(rel_vel, 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])
+            new_vels = _calc_post_col_vels(velocities[p1], velocities[p2], mass, mass, rel_vel, R[1], R[2])
             velocities[p1] = new_vels[0]
             velocities[p2] = new_vels[1]
-        
-        i += 2
-    
+
     return velocities
- 
-@njit   
+
+@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]:
+        if not number_children[i] and number_elements[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()
@@ -102,53 +102,54 @@ def clear(self):
         self.domain = None
         self.sigma_T = 3.631681e-19
         self.mass = None
-        
-    def advance(self, dt):
+
+    def advance(self, dt, collisions=True, octree=True):
         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)
+
+        if octree:
+            self.octree.build(self.particles.Pos)
+        if collisions and octree:
+            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, u = None):
         box = np.array(box)
         N = int(round(oc.get_V(box) * n / self.w))
         print("creating {} particles".format(N))
         self.particles.create_particles(box, self.mass, T, N)
-        
+
         if u is not None:
             velocities = self.particles.Vel
             u = np.array(u)
-            
+
             for i in range(len(velocities)):
                 velocities[i] += u
-            
+
             self.particles.VelPos = (velocities, self.particles.Pos)
-        
+
         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 f1dd8ad..18afcd4 100644
--- a/dsmc/octree.py
+++ b/dsmc/octree.py
@@ -1,6 +1,7 @@
 import numpy as np
 import numpy.typing as npt
 from numba import njit
+import numba as nb
 
 fmin = np.finfo(float).min
 fmax = np.finfo(float).max
@@ -128,6 +129,58 @@ def _create_boxes(box):
     child_geo8 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (box[2][0], half[2])))
     
     return [child_geo1, child_geo2, child_geo3, child_geo4, child_geo5, child_geo6, child_geo7, child_geo8]
+
+@njit
+def _get_min_aspect_ratio(box, axis):
+    half = np.array([0.5*(box[i][1] - box[i][0]) for i in range(3)])
+    
+    match axis:
+        case 0:
+            return min(half[0] / half[1], half[0] / half[2]);
+        case 1:
+            return min(half[1] / half[0], half[1] / half[2]);
+        case 2:
+            return min(half[2] / half[1], half[2] / half[0]);
+
+@njit
+def _devide(box, axis):
+    half = 0.5*(box[axis][0] + box[axis][1])
+    box1 = np.copy(box)
+    box2 = np.copy(box)
+    
+    box1[axis][0] = box[axis][0]
+    box1[axis][1] = half
+    
+    box2[axis][0] = half
+    box2[axis][1] = box[axis][1]
+    
+    return (box1, box2)
+
+@njit
+def _create_combined_boxes(box, min_aspect_ratio):
+    boxes = np.empty((15, 3, 2))
+    boxes[0] = box
+    N = 0
+    Nold = 0
+    q = 1
+    
+    for i in range(3):
+        if _get_min_aspect_ratio(box, i) > min_aspect_ratio:
+            for b in range(Nold, Nold + 2**N):
+                new_boxes = _devide(boxes[b], i)
+                boxes[q] = new_boxes[0]
+                boxes[q + 1] = new_boxes[1]
+                q += 2
+            Nold += 2**N
+            N += 1
+            
+    N = 2**N
+    new_boxes = np.empty((N, 3, 2))
+    
+    for b in range(N):
+        new_boxes[b] = boxes[Nold + b]
+            
+    return new_boxes
     
 class Leaf:
     def __init__(self):
@@ -141,6 +194,7 @@ def __init__(self):
 class Octree:
     def __init__(self):
         self.clear()
+        self.min_aspect_ratio = 2.0/3.0
         
     def clear(self):
         self.cell_boxes = []
@@ -179,10 +233,17 @@ def _create_root(self, positions):
         
     def _progress(self, leaf_id, positions):
         if _is_resolved(self.cell_boxes[leaf_id], self.leafs[leaf_id].number_elements, self.w, self.sigma_T, self.Nmin, self.Nmax):
-            self.leafs[leaf_id].number_children = 8
+            
             self.leafs[leaf_id].id_first_child = self.cell_offsets[-1]
-            self.cell_offsets[-1] += 8
-            self.cell_boxes += _create_boxes(self.cell_boxes[leaf_id])
+            
+            new_boxes = _create_combined_boxes(self.cell_boxes[leaf_id], self.min_aspect_ratio)
+            self.cell_offsets[-1] += len(new_boxes)
+            self.leafs[leaf_id].number_children = len(new_boxes)
+            
+            for box in new_boxes:
+                self.cell_boxes.append(box)
+            
+            #raise Exception()
         else:
             pass
            
diff --git a/dsmc/particles.py b/dsmc/particles.py
index 54723fc..7e709c3 100644
--- a/dsmc/particles.py
+++ b/dsmc/particles.py
@@ -37,7 +37,7 @@ def x2velocity(x, mass):
     -------
     speed of particle : float
     """
-    return math.sqrt((2.0/3.0) * x * kb /mass)
+    return math.sqrt(2.0 * x * kb /mass)
 
 @njit
 def get_vel(T, mass):
@@ -53,7 +53,8 @@ def get_vel(T, mass):
     -------
     velocity : np.array, shape = (3, 1)
     """
-    return np.array([(-1)**(int(2*np.random.random())) * x2velocity(box_muller(T), mass) for _ in range(3)])
+    v = np.random.random(3)
+    return v * x2velocity(box_muller(T), mass) / np.linalg.norm(v)
 
 @njit
 def get_velocities(T, mass, N):
diff --git a/dsmc/writer.py b/dsmc/writer.py
index 19a2a45..c181d62 100644
--- a/dsmc/writer.py
+++ b/dsmc/writer.py
@@ -1,7 +1,5 @@
-from . import octree as oc
-
-def write_buttom_leafs(octree):
-	f = open("octree.vtu", "w")
+def write_buttom_leafs(octree, file_name="octree.vtu"):
+	f = open(file_name, "w")
 	
 	_write_header(f)
 	_wrtie_body(f, octree)
diff --git a/examples/T_sample.py b/examples/T_sample.py
new file mode 100644
index 0000000..f40cb06
--- /dev/null
+++ b/examples/T_sample.py
@@ -0,0 +1,37 @@
+import dsmc
+import dsmc.particles as prt
+import matplotlib.pyplot as plt
+import numpy as np
+import math
+
+def maxwell(x, T):
+	return 2.0 * np.sqrt(x) * np.exp(-x/T) / (math.pow(T, 3.0/2.0) * np.sqrt(math.pi))
+    
+def calc_x(velocities, mass):
+	return np.array([mass*vel.dot(vel)/(2.0*prt.kb) for vel in velocities])
+    
+if __name__ == '__main__':
+	solver = dsmc.DSMC()
+	domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3))
+	w = 1e6
+	mass = 6.6422e-26
+	T = 300
+	n = 1e+20
+	Nbins = 100
+    
+	solver.set_domain(domain)
+	solver.set_weight(w)
+	solver.set_mass(mass)
+    
+	solver.create_particles(domain, T, n)
+    
+	print(prt.calc_temperature(solver.particles.Vel, mass))
+	
+	x = calc_x(solver.particles.Vel, solver.mass)
+
+	xm = np.linspace(0, 3500, 1000)
+	dist = [maxwell(xi, 300) for xi in xm]
+    
+	plt.plot(xm, dist)
+	plt.hist(x, Nbins, density=True)
+	plt.show()
diff --git a/examples/push_bound_test.py b/examples/push_bound_test.py
new file mode 100644
index 0000000..0c0798e
--- /dev/null
+++ b/examples/push_bound_test.py
@@ -0,0 +1,40 @@
+import dsmc
+import dsmc.octree as oc
+import dsmc.writer as wrt
+import numpy as np
+
+if __name__ == '__main__':
+	# general parameters
+    solver = dsmc.DSMC()
+    tree = oc.Octree()
+    domain = [(-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)]
+    positions = np.array([(-0.5, -0.5, -0.5), (0.5, 0.5, 0.5)])
+    dt = 1e-5
+    w = 1e+17
+    n = 1e18
+    mass = 6.6422e-26
+    T = 300
+    niter = 300
+    
+    # setup solver
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+    
+    solver.create_particles(domain, T, n)
+    
+    tree.build(positions)
+    
+    for it in range(niter):
+        print("iteration {:4}/{:4}".format(it + 1, niter), end="\r", flush=True)
+        solver.advance(dt, collisions=False, octree=False)
+        
+        with open("particles_{}.csv".format(it), "w") as file:
+            file.write("x, y, z\n")
+            for pos in solver.particles.Pos:
+                file.write("{},{},{}\n".format(pos[0], pos[1], pos[2]))
+                
+    wrt.write_buttom_leafs(tree, "octree.vtu")
+            
+    
+    print("done")
\ No newline at end of file
diff --git a/examples/shock_tube.py b/examples/shock_tube.py
index b10df19..37ef606 100644
--- a/examples/shock_tube.py
+++ b/examples/shock_tube.py
@@ -1,10 +1,8 @@
 import dsmc
 import dsmc.diagnostics as dia
-import matplotlib.pyplot as plt
 
-def write2file(solver, n_file, T_file, p_file):
+def write2file(solver, n_file, T_file, p_file, Nbins):
     bins, box, x = dia.sort_bin(solver.particles.Pos, 2, Nbins)
-    N = [len(b) for b in bins]
     n = dia.calc_n(bins, box, 2, solver.w)
     T = dia.calc_T(bins, solver.particles.Vel, mass)
     p = dia.calc_p(n, T)
@@ -54,7 +52,7 @@ def write2file(solver, n_file, T_file, p_file):
     for it in range(niter):
         print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
         solver.advance(dt)        
-        write2file(solver, n_file, T_file, p_file)
+        write2file(solver, n_file, T_file, p_file, Nbins)
 
 	# close files
     n_file.close()
diff --git a/examples/temp_relax.py b/examples/temp_relax.py
index bed76ac..dd369fb 100644
--- a/examples/temp_relax.py
+++ b/examples/temp_relax.py
@@ -3,6 +3,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import math
+import time
 
 def maxwell(x, T):
     return 2.0 * np.sqrt(x) * np.exp(-x/T) / (math.pow(T, 3.0/2.0) * np.sqrt(math.pi))
@@ -17,7 +18,7 @@ def calc_x(velocities, mass):
     dt = 1e-5
     w = 2.4134e+7
     mass = 6.6422e-26
-    niter = 500
+    niter = 100
     Nbins = 100
     
     # particles
@@ -33,17 +34,24 @@ def calc_x(velocities, mass):
     
     solver.create_particles(Box, T, n, u)
     
+    # time
+    start_time = time.time()
+    
     for it in range(niter):
         print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
         solver.advance(dt)
         x = calc_x(solver.particles.Vel, solver.mass)
 
-    xm = np.linspace(0, 3500, 1000)
-    dist = [maxwell(xi, 300) for xi in xm]
+    Tnew = prt.calc_temperature(solver.particles.Vel, solver.mass)
+    xm = np.linspace(0, 30000, 1000)
+    dist = [maxwell(xi, Tnew) for xi in xm]
+    
+    print("")
+    print("--- %s seconds ---" % (time.time() - start_time))
     
     plt.plot(xm, dist)
     plt.hist(x, Nbins, density=True)
-    plt.show()  
+    plt.title("T = {:.3f}K".format(Tnew))
+    plt.show() 
     
-    print("")
     print('done')
diff --git a/setup.py b/setup.py
index 005ff89..689ad62 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import setup
 setup(
     name='dsmc',
-    version='0.5.1',
+    version='0.6.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 90cdf10..bd7a9e2 100644
--- a/tests/unit/test_dsmc/dsmc.py
+++ b/tests/unit/test_dsmc/dsmc.py
@@ -9,13 +9,14 @@ def test_Constructor(self):
     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])
+        vel_rel = np.linalg.norm(vel1 - vel2)
         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)
+        res = ds._calc_prob(vel_rel, sigma_T, Vc, dt, w, N)
         
         self.assertEqual(np.linalg.norm(vel1 - vel2), res)
         
diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py
index 56a1874..5eb9a57 100644
--- a/tests/unit/test_dsmc/octree.py
+++ b/tests/unit/test_dsmc/octree.py
@@ -220,7 +220,103 @@ def test__create_boxes(self):
             V += oc.get_V(box)
             
         self.assertEqual(oc.get_V(box_orig), V)
+        
+    def test__get_min_aspect_ratio_1(self):
+        box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)])
+        axis1 = 0
+        axis2 = 1
+        axis3 = 2
+        
+        self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1))
+        self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2))
+        self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3))
+        
+    def test__get_min_aspect_ratio_2(self):
+        box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)])
+        axis1 = 0
+        axis2 = 1
+        axis3 = 2
+            
+        self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1))
+        self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2))
+        self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3))
+        
+    def test__devide(self):
+        box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)])
+        box_x1, box_x2 = oc._devide(box, 0)
+        box_y1, box_y2 = oc._devide(box, 1)
+        box_z1, box_z2 = oc._devide(box, 2)
+        
+        boxes = ((box_x1, box_x2), (box_y1, box_y2), (box_z1, box_z2))
+        
+        half = np.array([0.5*(box[i][0] + box[i][1]) for i in range(3)])
+        
+        for b in range(3):
+            box_a = boxes[b][0]
+            box_b = boxes[b][1]
+            
+            self.assertEqual(box_a[b][0], box[b][0])
+            self.assertEqual(box_a[b][1], half[b])
+            
+            self.assertEqual(box_b[b][0], half[b])
+            self.assertEqual(box_b[b][1], box[b][1])
+            
+            for i in range(3):
+                for j in range(2):
+                    if b != i:
+                        self.assertEqual(box_a[i][j], box[i][j])
+                        self.assertEqual(box_b[i][j], box[i][j])
+                        
+    def test__create_combined_boxes_1(self):
+        box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)])
+        min_aspect_ratio = 0.0
+        boxes_old = oc._create_boxes(box)
+        boxes_new = oc._create_combined_boxes(box, min_aspect_ratio)
+        V = 0.0
+        
+        for b in boxes_new:
+            V += oc.get_V(b)
             
+        self.assertEqual(oc.get_V(box), V)
+        self.assertEqual(len(boxes_old), len(boxes_new))
+        
+    def test__create_combined_boxes_2(self):
+        box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)])
+        min_aspect_ratio = 0.05
+        boxes_new = oc._create_combined_boxes(box, min_aspect_ratio)
+        V = 0.0
+        
+        for b in boxes_new:
+            V += oc.get_V(b)
+            
+        self.assertEqual(oc.get_V(box), V)
+        self.assertEqual(4, len(boxes_new))
+        
+    def test__create_combined_boxes_3(self):
+        box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)])
+        min_aspect_ratio = 0.05
+        boxes_new = oc._create_combined_boxes(box, min_aspect_ratio)
+        V = 0.0
+        
+        for b in boxes_new:
+            V += oc.get_V(b)
+            
+        self.assertEqual(oc.get_V(box), V)
+        self.assertEqual(4, len(boxes_new))
+        
+    def test__create_combined_boxes_4(self):
+        box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)])
+        min_aspect_ratio = 0.0
+        boxes_new = oc._create_combined_boxes(box, min_aspect_ratio)
+        V = 0.0
+        
+        for b in boxes_new:
+            V += oc.get_V(b)
+            
+        self.assertEqual(oc.get_V(box), V)
+        self.assertEqual(8, len(boxes_new))
+        
+        print(boxes_new)
             
 class TestOctreeOctree(unittest.TestCase):
     def test_build(self):