diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp
index 97bc986d..3ffa8b23 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.cpp
+++ b/image_analysis/geom_refinement/XtalOptimizer.cpp
@@ -9,7 +9,8 @@ struct XtalResidual {
     XtalResidual(double x, double y,
                  double lambda,
                  double pixel_size,
-                 const Eigen::Matrix3d &gonio_back_rot,
+                 Coord rot_axis,
+                 double angle_rad,
                  double exp_h, double exp_k,
                  double exp_l,
                  gemmi::CrystalSystem symmetry)
@@ -19,7 +20,8 @@ struct XtalResidual {
           exp_h(exp_h),
           exp_k(exp_k),
           exp_l(exp_l),
-          gonio_back_rot_(gonio_back_rot),
+          rot_axis(rot_axis),
+          angle_rad(angle_rad),
           symmetry(symmetry) {
     }
 
@@ -59,7 +61,15 @@ struct XtalResidual {
 
         // Apply goniometer "back-to-start" rotation:
         // brings observed reciprocal from image orientation into reference crystal frame
-        Eigen::Matrix<T, 3, 3> R_gonio_back = gonio_back_rot_.cast<T>();
+        T a[3];
+        a[0] = T(angle_rad*rot_axis.x);
+        a[1] = T(angle_rad*rot_axis.y);
+        a[2] = T(angle_rad*rot_axis.z);
+
+        T rot_arr[9];
+        ceres::AngleAxisToRotationMatrix(a, rot_arr);
+
+        Eigen::Matrix<T, 3, 3> R_gonio_back(rot_arr);
         Eigen::Matrix<T, 3, 1> e_obs_recip = R_gonio_back * e_obs_recip_raw;
 
         Eigen::Matrix<T, 3, 1> e_pred;
@@ -139,7 +149,8 @@ struct XtalResidual {
     const double exp_h;
     const double exp_k;
     const double exp_l;
-    const Eigen::Matrix3d gonio_back_rot_;
+    const Coord rot_axis;
+    const double angle_rad;
     gemmi::CrystalSystem symmetry;
 };
 
@@ -419,9 +430,12 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         for (const auto &pt: spots) {
             if (!data.index_ice_rings && pt.ice_ring)
                 continue;
-            
+
             Eigen::Matrix3d gonio_back_rot = Eigen::Matrix3d::Identity();
 
+            Coord axis = Coord(1,0,0);
+            float angle_rad = 0.0;
+
             Coord recip = data.geom.DetectorToRecip(pt.x, pt.y);
 
             if (data.axis) {
@@ -429,13 +443,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                 auto rot = data.axis->GetTransformationAngle(angle_deg);
                 recip = rot * recip;
 
-                Eigen::Matrix3d Rg;
-                for (int i = 0; i < 3; ++i) {
-                    for (int j = 0; j < 3; ++j) {
-                        Rg(i, j) = static_cast<double>(rot.arr()[i + j * 3]);
-                    }
-                }
-                gonio_back_rot = Rg.transpose();
+                angle_rad = angle_deg * M_PI / 180.0;
+                axis = data.axis->GetAxis();
             }
 
             double h_fp = recip * vec0;
@@ -456,7 +465,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                     new XtalResidual(pt.x, pt.y,
                                      data.geom.GetWavelength_A(),
                                      data.geom.GetPixelSize_mm(),
-                                     gonio_back_rot,
+                                     axis,
+                                     angle_rad,
                                      h, k, l,
                                      data.crystal_system)),
                 nullptr,
@@ -547,7 +557,9 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
             data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
 
         if (!data.refine_unit_cell) {
-            Coord rot_vector{static_cast<float>(latt_vec0[0]), static_cast<float>(latt_vec0[1]), static_cast<float>(latt_vec0[2])};
+            Coord rot_vector{
+                static_cast<float>(latt_vec0[0]), static_cast<float>(latt_vec0[1]), static_cast<float>(latt_vec0[2])
+            };
             Coord rot_vector_norm = rot_vector.Normalize();
             std::cout << "X " << rot_vector.Length() * 180.0 / M_PI << " " << rot_vector_norm << std::endl;
         }