ScaleOnTheFly: Calculate rotation partiality correctly

image_analysis/scale_merge/ScaleOnTheFly.cpp

@@ -13,6 +13,17 @@ namespace {
         return 1.0 / x;
     }
 
+    float RotationPartiality( double delta_phi_deg,
+                     double zeta,
+                     double mosaicity_deg,
+                     double wedge_deg) {
+        const double half_wedge = wedge_deg / 2.0;
+        const double c1 = zeta / std::sqrt(2.0);
+        const double arg_plus = (delta_phi_deg + half_wedge) * c1 / mosaicity_deg;
+        const double arg_minus = (delta_phi_deg - half_wedge) * c1 / mosaicity_deg;
+        return static_cast<float>((std::erf(arg_plus) - std::erf(arg_minus)) / 2.0);
+    }
+
     class ScalingResidual {
     protected:
         const double Iobs;
@@ -204,10 +215,8 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &reflections, s
                 r.partiality = 1.0;
                 break;
             case PartialityModel::Rotation: {
-                double partiality = 0.0;
-                ScalingRotationResidual res(r, 0, 0);
-                if (res(&result.G, &result.B, &result.mos, &result.wedge, &partiality))
-                    r.partiality = static_cast<float>(partiality);
+                if (std::isfinite(r.delta_phi_deg) && std::isfinite(r.zeta) && result.mos > 1e-6)
+                    r.partiality = RotationPartiality(r.delta_phi_deg, r.zeta, result.mos, result.wedge);
                 break;
             }
             default: