Skip to content
Snippets Groups Projects
Commit a30aca5a authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

[386a] Repair slicing and ff computation for HorizontalCylinder (#386) ()

Merging branch '386a'  into 'main'.

See merge request !1228
parents 52f7a433 c622a88b
No related branches found
No related tags found
1 merge request!1228Repair slicing and ff computation for HorizontalCylinder (#386)
Pipeline #84632 passed
......@@ -75,8 +75,8 @@ IFormFactor* doSlice(const IFormFactor* ff, double dz_bottom, double dz_top)
f->height() - dz_bottom - dz_top);
}
if (const auto* f = dynamic_cast<const HorizontalCylinder*>(ff)) {
return new HorizontalCylinder(f->radius(), f->length(), -f->radius() + dz_bottom,
f->radius() - dz_top);
return new HorizontalCylinder(f->radius(), f->length(), f->slice_bottom() + dz_bottom,
f->slice_top() - dz_top);
}
if (const auto* f = dynamic_cast<const Sphere*>(ff)) {
return new TruncatedSphere(f->radius(), f->height() - dz_bottom, dz_top);
......
......@@ -57,14 +57,14 @@ complex_t HorizontalCylinder::formfactor(C3 q) const
radial_part = M_TWOPI * R * R * Math::Bessel::J1c(qR) * exp_I(q.z() * R);
else
radial_part =
2.
* ComplexIntegrator().integrate(
[=](double z) {
double y = sqrt(R * R - z * z);
return y * Math::sinc(q.y() * y) * exp_I(q.z() * (z - m_slice_bottom));
},
m_slice_bottom, m_slice_top);
// integration variable substituted as z = R * sin(t)
radial_part = 2. * pow(R, 2) * exp_I(q.z() * (-m_slice_bottom))
* ComplexIntegrator().integrate(
[=](double t) {
return pow(cos(t), 2) * Math::sinc(q.y() * R * cos(t))
* exp_I(q.z() * R * sin(t));
},
asin(m_slice_bottom / R), asin(m_slice_top / R));
return radial_part * axial_part;
}
......
......@@ -41,20 +41,20 @@ def get_sample(particle, sliced):
class SlicedSpheresTest(unittest.TestCase):
def runPlainFF(self, ff):
def runPlainFF(self, ff, eps = 1e-14):
particle = ba.Particle(matParticle, ff)
diff = infrastruct.diff_MiniGISAS(get_sample(particle, False),
get_sample(particle, True))
self.assertLess(diff, 1e-9)
self.assertLess(diff, eps)
def testSlicingPlainFF(self):
self.runPlainFF(ba.Cone(8., 9., 80*ba.deg))
self.runPlainFF(ba.Cylinder(3., 9.))
self.runPlainFF(ba.EllipsoidalCylinder(3., 4., 9.))
# self.runPlainFF(ba.HemiEllipsoid(7., 8., 9.)) # yet unsupported
self.runPlainFF(ba.HorizontalCylinder(5., 19.)) # TODO restore tol=1e-13
#self.runPlainFF(ba.HorizontalCylinder(5., 19., -4., 4.)) TODO wrong!
#self.runPlainFF(ba.HorizontalCylinder(6., 19., -3., 3.))
self.runPlainFF(ba.HorizontalCylinder(5., 19.), 5e-12)
self.runPlainFF(ba.HorizontalCylinder(5., 19., -4., 4.), 5e-12)
self.runPlainFF(ba.HorizontalCylinder(6., 19., -3., 3.), 5e-12)
self.runPlainFF(ba.Sphere(5.))
self.runPlainFF(ba.Spheroid(4., 4.))
self.runPlainFF(ba.TruncatedSphere(7., 9., 1.))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment