diff --git a/src/region_ellipsoid.cpp b/src/region_ellipsoid.cpp
index a413f584fc..24e9386f61 100644
--- a/src/region_ellipsoid.cpp
+++ b/src/region_ellipsoid.cpp
@@ -198,36 +198,30 @@ int RegEllipsoid::surface_interior(double *x, double cutoff)
 
     if (r_r > rc_r) {
       // sort the values
-      int sorting[3] = {0, 1, 2};
-      double axes[3] = {c, b, a};
-      double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
+      double axes[3] = {a,b,c};
+      double coords[3] = {fabs(x[0] - xc), fabs(x[1] - yc), fabs(x[2] - zc)};
 
-      if (axes[1] < axes[0]) {
-        sorting[0] = 1;
-        sorting[1] = 0;
-        axes[0] = b;
-        axes[1] = c;
-      }
-
-      if (axes[2] < axes[1]) {
-        int ti = sorting[2];
-        sorting[2] = sorting[1];
-        sorting[1] = ti;
-        double td = axes[2];
-        axes[2] = axes[1];
-        axes[1] = td;
-        if (axes[1] < axes[0]) ti = sorting[1];
-        sorting[1] = sorting[0];
-        sorting[0] = ti;
-        td = axes[1];
-        axes[1] = axes[0];
-        axes[0] = td;
-      }
+     int min, max;
+     if (axes[0] < axes[1]) {
+       min = 0;
+       max = 1;
+     } else {
+       min = 1;
+       max = 0;
+     }
+     if (axes[min] > axes[2]) {
+       min = 2;
+     }
+     if (axes[max] < axes[2]) {
+       max = 2;
+     }
+     int mid = 3 - min - max;
+     int sorting[3] = {min, mid, max};
 
       double x0[3];
       contact[0].r =
-          DistancePointEllipsoid(axes[2], axes[1], axes[0], coords[sorting[2]], coords[sorting[1]],
-                                 coords[sorting[0]], x0[2], x0[1], x0[0]);
+          DistancePointEllipsoid(axes[sorting[2]], axes[sorting[1]], axes[sorting[0]], coords[sorting[2]], coords[sorting[1]],
+                                 coords[sorting[0]], x0[sorting[2]], x0[sorting[1]], x0[sorting[0]]);
       contact[0].delx = x[0] - (copysign(x0[sorting[2]], x[0] - xc) + xc);
       contact[0].dely = x[1] - (copysign(x0[sorting[1]], x[1] - yc) + yc);
       contact[0].delz = x[2] - (copysign(x0[sorting[0]], x[2] - zc) + zc);
@@ -299,36 +293,30 @@ int RegEllipsoid::surface_exterior(double *x, double cutoff)
 
     if (r_r < rc_r) {
       // sort the values
-      int sorting[3] = {0, 1, 2};
-      double axes[3] = {c, b, a};
-      double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
+      double axes[3] = {a,b,c};
+      double coords[3] = {fabs(x[0] - xc), fabs(x[1] - yc), fabs(x[2] - zc)};
 
-      if (axes[1] < axes[0]) {
-        sorting[0] = 1;
-        sorting[1] = 0;
-        axes[0] = b;
-        axes[1] = c;
+      int min, max;
+      if (axes[0] < axes[1]) {
+        min = 0;
+        max = 1;
+      } else {
+        min = 1;
+        max = 0;
       }
-
-      if (axes[2] < axes[1]) {
-        int ti = sorting[2];
-        sorting[2] = sorting[1];
-        sorting[1] = ti;
-        double td = axes[2];
-        axes[2] = axes[1];
-        axes[1] = td;
-        if (axes[1] < axes[0]) ti = sorting[1];
-        sorting[1] = sorting[0];
-        sorting[0] = ti;
-        td = axes[1];
-        axes[1] = axes[0];
-        axes[0] = td;
+      if (axes[min] > axes[2]) {
+        min = 2;
       }
+      if (axes[max] < axes[2]) {
+        max = 2;
+      }
+      int mid = 3 - min - max;
+      int sorting[3] = {min, mid, max};
 
       double x0[3];
       contact[0].r =
-          DistancePointEllipsoid(axes[2], axes[1], axes[0], coords[sorting[2]], coords[sorting[1]],
-                                 coords[sorting[0]], x0[2], x0[1], x0[0]);
+          DistancePointEllipsoid(axes[sorting[2]], axes[sorting[1]], axes[sorting[0]], coords[sorting[2]], coords[sorting[1]],
+                                 coords[sorting[0]], x0[sorting[2]], x0[sorting[1]], x0[sorting[0]]);
       contact[0].delx = x[0] - (copysign(x0[sorting[2]], x[0] - xc) + xc);
       contact[0].dely = x[1] - (copysign(x0[sorting[1]], x[1] - yc) + yc);
       contact[0].delz = x[2] - (copysign(x0[sorting[0]], x[2] - zc) + zc);