diff --git a/src/kanapy/api.py b/src/kanapy/api.py
index 6e77b7e1..d245b56c 100644
--- a/src/kanapy/api.py
+++ b/src/kanapy/api.py
@@ -372,7 +372,11 @@ def plot_ellipsoids(self, cmap='prism', dual_phase=False):
         """ Generates plot of particles"""
         if self.particles is None:
             raise ValueError('No particle to plot. Run pack first.')
-        plot_ellipsoids_3D(self.particles, cmap=cmap, dual_phase=dual_phase)
+        hmin = min(self.rve.size)
+        asp_arr = [int(self.rve.size[0] / hmin),
+                   int(self.rve.size[1] / hmin),
+                   int(self.rve.size[2] / hmin)]
+        plot_ellipsoids_3D(self.particles, cmap=cmap, dual_phase=dual_phase, asp_arr=asp_arr)
 
     def plot_particles(self, cmap='prism', dual_phase=False, plot_hull=True):
         """ Generates plot of particles"""
@@ -380,8 +384,12 @@ def plot_particles(self, cmap='prism', dual_phase=False, plot_hull=True):
             raise ValueError('No particle to plot. Run pack first.')
         if self.particles[0].inner is None:
             raise ValueError('Ellipsoids without inner polygon cannot be plotted.')
+        hmin = min(self.rve.size)
+        asp_arr = [int(self.rve.size[0] / hmin),
+                   int(self.rve.size[1] / hmin),
+                   int(self.rve.size[2] / hmin)]
         plot_particles_3D(self.particles, cmap=cmap,
-                          dual_phase=dual_phase, plot_hull=plot_hull)
+                          dual_phase=dual_phase, plot_hull=plot_hull, asp_arr=asp_arr)
 
     def plot_voxels(self, sliced=False, dual_phase=False, cmap='prism', ori=None,
                     color_key=0, silent=False):
@@ -420,9 +428,13 @@ def plot_voxels(self, sliced=False, dual_phase=False, cmap='prism', ori=None,
                 clist = None
         else:
             clist = None
-
+        hmin = min(self.rve.size)
+        asp_arr = [int(self.rve.size[0] / hmin),
+                   int(self.rve.size[1] / hmin),
+                   int(self.rve.size[2] / hmin)]
         fig = plot_voxels_3D(data, Ngr=np.sum(self.ngrains), sliced=sliced,
-                             dual_phase=dual_phase, cmap=cmap, clist=clist, silent=silent)
+                             dual_phase=dual_phase, cmap=cmap, clist=clist,
+                             silent=silent, asp_arr=asp_arr)
         if silent:
             return fig
 
@@ -433,8 +445,12 @@ def plot_grains(self, geometry=None, cmap='prism', alpha=0.4,
             geometry = self.geometry
         if geometry is None:
             raise ValueError('No polygons for grains defined. Run generate_grains() first')
+        hmin = min(self.rve.size)
+        asp_arr = [int(self.rve.size[0] / hmin),
+                   int(self.rve.size[1] / hmin),
+                   int(self.rve.size[2] / hmin)]
         plot_polygons_3D(geometry, cmap=cmap, alpha=alpha, ec=ec,
-                         dual_phase=dual_phase)
+                         dual_phase=dual_phase, asp_arr=asp_arr)
 
     def plot_stats(self, data=None,
                    gs_data=None, gs_param=None,
diff --git a/src/kanapy/initializations.py b/src/kanapy/initializations.py
index 3f2ec8e2..2c92565b 100644
--- a/src/kanapy/initializations.py
+++ b/src/kanapy/initializations.py
@@ -101,7 +101,7 @@ def __init__(self, stats_dicts, nsteps=1000, from_voxels=False, poly=None):
             Number of alloy in ICAMS CP-UMAT
         """
 
-        def init_particles():
+        def init_particles(ip):
             """
             Extract statistical microstructure information from data dictionary
             and initalize particles for packing accordingly
@@ -135,7 +135,8 @@ def gen_data_basic(pdict):
                 K = individualK / np.sum(individualK)
 
                 # Total number of ellipsoids for packing density given by volume fraction
-                num_f = np.divide(K * phase_vf[-1] * np.prod(self.size), volume_array)
+                ip = pdict['Phase']
+                num_f = np.divide(K * phase_vf[ip] * np.prod(self.size), volume_array)
                 """print(f'Particle numbers: {num_f}')
                 print(f'Total volume of particles: {np.sum(np.multiply(num_f, volume_array))}')"""
                 # Rounding particle numbers to integer values keeping total volume constant
@@ -151,7 +152,7 @@ def gen_data_basic(pdict):
                     num[i] = v1
                 """print(f'Volume after rounding: {np.sum(np.multiply(num, volume_array))}')
                 print(f'Particle numbers: {num}')
-                print(f'Phase volume fraction: {phase_vf[-1]}')"""
+                print(f'Phase volume fraction: {phase_vf[ip]}')"""
                 totalEllipsoids = int(np.sum(num))
 
                 # Duplicate the diameter values
@@ -389,7 +390,7 @@ def gen_data_elong(pdict):
             if not from_voxels:
                 if stats["Grain type"] not in ["Elongated", "Equiaxed", "Free"]:
                     raise ValueError('The value for "Grain type" must be either "Equiaxed" or "Elongated".')
-                part_dict = init_particles()
+                part_dict = init_particles(ip)
                 self.particle_data.append(part_dict)
                 self.nparticles.append(part_dict['Number'])
         print('  RVE characteristics:')
diff --git a/src/kanapy/plotting.py b/src/kanapy/plotting.py
index a48e350b..56b4526a 100644
--- a/src/kanapy/plotting.py
+++ b/src/kanapy/plotting.py
@@ -37,7 +37,7 @@
 
 def plot_voxels_3D(data, Ngr=1, sliced=False, dual_phase=False,
                    mask=None, cmap='prism', alpha=1.0, silent=False,
-                   clist=None):
+                   clist=None, asp_arr=None):
     """
     Plot voxels in microstructure, each grain with a different color. Sliced
     indicates whether one eighth of the box should be removed to see internal
@@ -75,7 +75,8 @@ def plot_voxels_3D(data, Ngr=1, sliced=False, dual_phase=False,
     Nx = data.shape[0]
     Ny = data.shape[1]
     Nz = data.shape[2]
-
+    if asp_arr is None:
+        asp_arr = [1, 1, 1]
     if mask is None:
         vox_b = np.full(data.shape, True, dtype=bool)
     else:
@@ -116,6 +117,7 @@ def plot_voxels_3D(data, Ngr=1, sliced=False, dual_phase=False,
     ax.set_xlim(right=Nx)
     ax.set_ylim(top=Ny)
     ax.set_zlim(top=Nz)
+    ax.set_box_aspect(asp_arr)
     if silent:
         return fig
     else:
@@ -123,7 +125,7 @@ def plot_voxels_3D(data, Ngr=1, sliced=False, dual_phase=False,
 
 
 def plot_polygons_3D(geometry, cmap='prism', alpha=0.4, ec=[0.5, 0.5, 0.5, 0.1],
-                     dual_phase=False, silent=False):
+                     dual_phase=False, silent=False, asp_arr=None):
     """
     Plot triangularized convex hulls of grains, based on vertices, i.e.
     connection points of 4 up to 8 grains or the end points of triple or quadruple
@@ -148,6 +150,8 @@ def plot_polygons_3D(geometry, cmap='prism', alpha=0.4, ec=[0.5, 0.5, 0.5, 0.1],
     None.
 
     """
+    if asp_arr is None:
+        asp_arr = [1, 1, 1]
     grains = geometry['Grains']
     pts = geometry['Points']
     Ng = np.amax(list(grains.keys()))
@@ -171,6 +175,7 @@ def plot_polygons_3D(geometry, cmap='prism', alpha=0.4, ec=[0.5, 0.5, 0.5, 0.1],
     ax.set(xlabel='x', ylabel='y', zlabel='z')
     ax.set_title('Polygonized microstructure')
     ax.view_init(30, 30)
+    ax.set_box_aspect(asp_arr)
     if silent:
         return fig
     else:
@@ -178,7 +183,7 @@ def plot_polygons_3D(geometry, cmap='prism', alpha=0.4, ec=[0.5, 0.5, 0.5, 0.1],
 
 
 
-def plot_ellipsoids_3D(particles, cmap='prism', dual_phase=False, silent=False):
+def plot_ellipsoids_3D(particles, cmap='prism', dual_phase=False, silent=False, asp_arr=None):
     """
     Display ellipsoids after packing procedure
     Parameters
@@ -190,7 +195,8 @@ def plot_ellipsoids_3D(particles, cmap='prism', dual_phase=False, silent=False):
     dual_phase : bool, optional
         Whether to display the ellipsoids in red/green contrast or in colors
     """
-
+    if asp_arr is None:
+        asp_arr = [1, 1, 1]
     fig = plt.figure(figsize=(6, 6))
     ax = fig.add_subplot(111, projection='3d')
     ax.set(xlabel='x', ylabel='y', zlabel='z')
@@ -211,6 +217,7 @@ def plot_ellipsoids_3D(particles, cmap='prism', dual_phase=False, silent=False):
         y = (pts[:, 1] + ctr[1]).reshape((100, 100))
         z = (pts[:, 2] + ctr[2]).reshape((100, 100))
         ax.plot_surface(x, y, z, rstride=4, cstride=4, color=color, linewidth=0, antialiased=False)
+    ax.set_box_aspect(asp_arr)
     if silent:
         return fig
     else:
@@ -218,7 +225,8 @@ def plot_ellipsoids_3D(particles, cmap='prism', dual_phase=False, silent=False):
 
 
 
-def plot_particles_3D(particles, cmap='prism', dual_phase=False, plot_hull=True, silent=False):
+def plot_particles_3D(particles, cmap='prism', dual_phase=False, plot_hull=True,
+                      silent=False, asp_arr=None):
     """
     Display inner polyhedra of ellipsoids after packing procedure
 
@@ -238,6 +246,8 @@ def plot_particles_3D(particles, cmap='prism', dual_phase=False, plot_hull=True,
     None.
 
     """
+    if asp_arr is None:
+        asp_arr = [1, 1, 1]
     fig = plt.figure(figsize=plt.figaspect(1))
     ax = fig.add_subplot(111, projection='3d')
     ax.set(xlabel='x', ylabel='y', zlabel='z')
@@ -272,6 +282,7 @@ def plot_particles_3D(particles, cmap='prism', dual_phase=False, plot_hull=True,
             y = (pts[:, 1] + ctr[None, 1]).reshape((100, 100))
             z = (pts[:, 2] + ctr[None, 2]).reshape((100, 100))
             ax.plot_surface(x, y, z, rstride=4, cstride=4, color=col, linewidth=0)
+    ax.set_box_aspect(asp_arr)
     if silent:
         return fig
     else:
diff --git a/tests/test_packing.py b/tests/test_packing.py
index 43fadaa2..2cef705d 100644
--- a/tests/test_packing.py
+++ b/tests/test_packing.py
@@ -16,7 +16,7 @@ def par_sim(mocker):
                'Minor_diameter1': np.array([2.15, 3.6, 1.15]),
                'Minor_diameter2': np.array([2.15, 3.6, 1.15]),
                'Tilt angle': np.array([92, 89.3, 85]),
-               'Phase': [0, 0, 0]}
+               'Phase': 0}
 
     sb = mocker.MagicMock()
     # Define attributes to mocker object
@@ -119,7 +119,7 @@ def test_packingRoutine():
           'Minor_diameter1': [1.5, 1.5],
           'Minor_diameter2': [1.5, 1.5],
           'Tilt angle': [86, 92],
-          'Phase': [0, 0]}
+          'Phase': 0}
     sb = Simulation_Box((3, 3, 3))
 
     # Test if the 'particle_generator' function is called once