diff --git a/Sample/HardParticle/HemiEllipsoid.cpp b/Sample/HardParticle/HemiEllipsoid.cpp
index 9b28c74c756fa862c191ee98431ff77ad86ef228..640e90a487808393cf434d6e5670d23d1d07477c 100644
--- a/Sample/HardParticle/HemiEllipsoid.cpp
+++ b/Sample/HardParticle/HemiEllipsoid.cpp
@@ -46,7 +46,8 @@ complex_t HemiEllipsoid::formfactor_at_bottom(C3 q) const
     const double R = m_radius_x;
     const double W = m_radius_y;
     const double H = m_height;
-    const complex_t G = std::sqrt((q.x() * R)*(q.x() * R) +  (q.y() * W) * (q.y() * W));
+    const complex_t G = std::sqrt((q.x() * R) * (q.x() * R) + (q.y() * W) * (q.y() * W));
+    const complex_t HQ = q.z() * H;
 
     if (std::abs(q.mag()) <= std::numeric_limits<double>::epsilon())
         return M_TWOPI * R * W * H / 3.;
@@ -54,10 +55,8 @@ complex_t HemiEllipsoid::formfactor_at_bottom(C3 q) const
     return M_TWOPI * H * R * W
            * ComplexIntegrator().integrate(
                [=](double z) {
-                   complex_t gamma = G * std::sqrt(1.0 - z * z);
-                   complex_t J1_gamma_div_gamma = Math::Bessel::J1c(gamma);
-
-                   return (1. - z*z) * J1_gamma_div_gamma * exp_I(q.z() * H * z);
+                   const double cz = 1 - z * z;
+                   return cz * Math::Bessel::J1c(G * std::sqrt(cz)) * exp_I(HQ * z);
                },
                0., 1.);
 }