#include "Base/Const/Units.h"

#include "Base/Type/Span.h"
#include "Sample/HardParticle/HardParticles.h"
#include "Sample/Scattering/Rotations.h"
#include "Tests/GTestWrapper/google_test.h"
#include <numbers>

using std::numbers::pi;
using namespace std::complex_literals;

void test_eps_q(const IFormfactor* p, C3 qdir, double eps)
{
    C3 q = eps * qdir;
    complex_t ff = p->formfactor(q);
    // std::cout<<"q="<<q<<" -> "<<std::setprecision(16)<<" ff0="<<V<<", ff ="<<ff<<"\n";
    double V = p->volume();
    double R = p->radialExtension();
    EXPECT_LE(real(ff), V * (1 + 4e-16));
    EXPECT_GT(real(ff), V * (1 - std::max(3e-16, 2 * eps * R * eps * R)));
    EXPECT_LT(std::abs(imag(ff)), 2 * eps * V * R);
}

void test_small_q(const IFormfactor* p, complex_t x, complex_t y, complex_t z)
{
    C3 q(x, y, z);
    test_eps_q(p, q, 1e-14);
    test_eps_q(p, q, 1e-11);
    test_eps_q(p, q, 1e-8);
    test_eps_q(p, q, 1e-5);
}

void test_ff(const IFormfactor* p)
{
    complex_t ff0 = p->formfactor(C3(0., 0., 0.));
    EXPECT_EQ(imag(ff0), 0.);
    double V = real(ff0);
    EXPECT_NEAR(p->volume(), V, 3e-16 * V);

    double R = p->radialExtension();
    if (R * R * R < V / 20 || R * R * R > 20 * V) {
        std::cerr << "WARNING: very elongated particle, or wrong radius; some tests disabled\n";
        std::cerr << std::setprecision(16) << "  V  =" << V << "\n  R^3=" << R * R * R << "\n";
    } else {
        test_small_q(p, 1, 0, 0);
        test_small_q(p, -1, 0, 0);
        test_small_q(p, 0, 1, 0);
        test_small_q(p, 0, -1, 0);
        test_small_q(p, 0, 0, 1);
        test_small_q(p, 0, 0, -1);
        test_small_q(p, 1, 1, 0);
        test_small_q(p, 0, 1, 1);
        test_small_q(p, 1, 0, 1);
        test_small_q(p, 1, -1, 0);
        test_small_q(p, 0, 1, -1);
        test_small_q(p, -1, 0, 1);
        test_small_q(p, 1, 1, 1);
        test_small_q(p, .7, .8, .9);
    }
    IFormfactor* clone = p->clone();
    EXPECT_EQ(clone->volume(), V);
    C3 q(.1, .2, complex_t(.3, .004));
    EXPECT_EQ(clone->formfactor(q), p->formfactor(q));
    delete clone;
}

const C3 q0 = {0.23, 0.11, -0.33};
const C3 q1 = {-0.52, -0.68, -0.69};
const C3 q2 = {0.16, -0.48, -0.46};
const C3 q3 = {-0.27, 0.87, 0.41};
const C3 q4 = {0.31, 0.36, -0.43};
const C3 q5 = {-0.60, 0.79, -0.37};
const C3 q6 = {-0.20, -0.89, 0.27};
const C3 q7 = {-0.18, 0.47, -0.06};
const C3 q8 = {-0.56, -0.94, 0.12};
const C3 q9 = {-0.51, -0.16, -0.42};

void test_ff_print(const IFormfactor* p, const C3 a, const C3 b, const C3 c,
                   std::complex<double> exp_a, std::complex<double> exp_b,
                   std::complex<double> exp_c)
{
    std::complex<double> ff_a = p->formfactor(a);
    std::complex<double> ff_b = p->formfactor(b);
    std::complex<double> ff_c = p->formfactor(c);
    std::cout << std::fixed;
    std::cout << "Formfactor Calculation for " << p->className() << "\n";
    std::cout << std::setprecision(2) << "F(" << a << ") = " << std::setprecision(15) << ff_a
              << "\n";
    std::cout << std::setprecision(2) << "F(" << b << ") = " << std::setprecision(15) << ff_b
              << "\n";
    std::cout << std::setprecision(2) << "F(" << c << ") = " << std::setprecision(15) << ff_c
              << "\n\n";
    double abs_a = abs(exp_a);
    EXPECT_NEAR(fabs(ff_a.real() - exp_a.real()), 0, abs_a * 1E-13);
    EXPECT_NEAR(fabs(ff_a.imag() - exp_a.imag()), 0, abs_a * 1E-13);
    double abs_b = abs(exp_b);
    EXPECT_NEAR(fabs(ff_b.real() - exp_b.real()), 0, abs_b * 1E-13);
    EXPECT_NEAR(fabs(ff_b.imag() - exp_b.imag()), 0, abs_b * 1E-13);
    double abs_c = abs(exp_c);
    EXPECT_NEAR(fabs(ff_c.real() - exp_c.real()), 0, abs_c * 1E-13);
    EXPECT_NEAR(fabs(ff_c.imag() - exp_c.imag()), 0, abs_c * 1E-13);
}

TEST(FormfactorBasic, Name)
{
    Box ff(1., 2., 3.);
    EXPECT_EQ(ff.className(), "Box");
    EXPECT_EQ(ff.className(), "Box");
}

TEST(FormfactorBasic, Pyramid2)
{
    EXPECT_THROW(Pyramid2(12, 14, -5, .8), std::runtime_error);
    EXPECT_THROW(Pyramid2(0, 1, 1, pi / 4), std::runtime_error);
    EXPECT_THROW(Pyramid2(1, 1, 1, 3 * pi / 4), std::runtime_error);

    double length = 12.;
    double width = 14.;
    double height = 5.;
    double alpha = 0.8;
    double tga = std::tan(alpha);
    double volume = height
                    * (length * width - (length + width) * height / tga
                       + 4.0 / 3.0 * height * height / (tga * tga));

    Pyramid2 particle(length, width, height, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q5, q3, q8, {-3.326507178430593, -4.221250124559902},
                  {-28.875234172319271, -14.799465747269439},
                  {-8.434962366982653, -1.380171136944693});
}

TEST(FormfactorBasic, Bar)
{
    EXPECT_THROW(BarGauss(-100, -20, 4), std::runtime_error);
    EXPECT_THROW(BarGauss(100, 20, -4), std::runtime_error);
    // TODO EXPECT_THROW(BarGauss(0, 20, 4), std::runtime_error);
    // TODO EXPECT_THROW(BarGauss(100, 0, 4), std::runtime_error);
    // TODO EXPECT_THROW(BarGauss(100, 20, 0), std::runtime_error);
    // test_ff_print(&particle,q0,q4,q1,0,0,0);

    double width = 20.;
    double height = 4.;
    double length = 100.0;
    double volume = height * width * length;

    BarGauss particle(length, width, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);
}

TEST(FormfactorBasic, Box)
{
    EXPECT_THROW(Box(-6, 5, -7), std::runtime_error);
    EXPECT_THROW(Box(6, 0, 7), std::runtime_error);

    double length = 6.;
    double height = 5.;
    double width = 7.;
    double volume = length * height * width;

    Box particle(length, width, height);
    EXPECT_NEAR(sqrt(pow(length / 2, 2) + pow(height / 2, 2) + pow(width / 2, 2)),
                particle.radialExtension(), 1e-14);
    EXPECT_DOUBLE_EQ(volume, particle.volume());

    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    EXPECT_NEAR(0., particle.spanZ(new RotationZ(17 * Units::deg)).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(new RotationZ(39 * Units::deg)).hig(), height * 1e-15);

    EXPECT_NEAR(-width / 2, particle.spanZ(new RotationX(90 * Units::deg)).low(), 1e-12);
    EXPECT_NEAR(-length / 2, particle.spanZ(new RotationY(90 * Units::deg)).low(), 1e-12);

    test_ff(&particle);
    test_ff_print(&particle, q9, q7, q2, {53.471450007876584, -93.217597445730974},
                  {119.400684160766147, -18.045648436947012},
                  {38.756783058601833, -86.601913488017942});
}

TEST(FormfactorBasic, CantellatedCube)
{
    EXPECT_THROW(CantellatedCube(0, 0), std::runtime_error);
    EXPECT_THROW(CantellatedCube(10, -1), std::runtime_error);
    EXPECT_THROW(CantellatedCube(10, 6), std::runtime_error);

    double L = 10.;
    double t = 2.; // side length of removed trirectangular tetrahedron at each vertex
    double volume = L * L * L - 6 * L * t * t + 16 * t * t * t / 3;

    CantellatedCube particle(L, t);
    EXPECT_DOUBLE_EQ(particle.volume(), volume);
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), L * 1e-15);
    EXPECT_NEAR(L, particle.spanZ(&rot).hig(), L * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q5, q6, {-33.391098526973408, -420.703291952477457},
                  {10.070159899276099, 35.125317795931998},
                  {-20.124862756519551, -89.660726461058985});
}

TEST(FormfactorBasic, Cone)
{
    EXPECT_THROW(Cone(-6, -5, 0.8), std::runtime_error);
    // TODO EXPECT_THROW(Cone(6, 5, 0), std::runtime_error);
    EXPECT_THROW(Cone(6, 5, 2 * pi / 3), std::runtime_error);

    double radius = 6.;
    double height = 5.;
    double alpha = 0.8;
    double tga = std::tan(alpha);
    double HdivRtga = height / tga / radius;
    double volume = pi / 3. * tga * radius * radius * radius
                    * (1. - (1. - HdivRtga) * (1. - HdivRtga) * (1. - HdivRtga));

    Cone particle(radius, height, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q8, q4, {154.836630212175464, -89.616777161335065},
                  {-7.158085433517508, 2.913985277509306},
                  {71.266836187605691, -75.316792284104494});
}

TEST(FormfactorBasic, Pyramid6)
{
    EXPECT_THROW(Pyramid6(-6, -5, 0.8), std::runtime_error);
    // TODO EXPECT_THROW(Pyramid6(6, 5, 0), std::runtime_error);
    EXPECT_THROW(Pyramid6(6, 5, 2 * pi / 3), std::runtime_error);

    double base_edge = 6.;
    double height = 5.;
    double alpha = 0.8;
    double tga = std::tan(alpha);
    double HdivRtga = 2. * height / tga / base_edge / std::sqrt(3.);
    double volume = 3. / 4. * tga * base_edge * base_edge * base_edge
                    * (1. - (1. - HdivRtga) * (1. - HdivRtga) * (1. - HdivRtga));

    Pyramid6 particle(base_edge, height, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());

    test_ff(&particle);
    test_ff_print(&particle, q9, q3, q5, {56.580851533712433, -49.055032503649102},
                  {-0.620299369844318, 17.442172408792921},
                  {-5.221867418516533, -12.118764972593246});
}

TEST(FormfactorBasic, CosineRippleBox)
{
    EXPECT_THROW(CosineRippleBox(-20, 4, -100), std::runtime_error);
    // TODO EXPECT_THROW(CosineRippleBox(0, 4, 100), std::runtime_error);
    // TODO EXPECT_THROW(CosineRippleBox(20, 0, 100), std::runtime_error);
    // TODO EXPECT_THROW(CosineRippleBox(20, 4, 0), std::runtime_error);

    double width = 20.;
    double height = 4.;
    double length = 100.0;
    double volume = 0.5 * height * width * length;

    CosineRippleBox particle(length, width, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q2, q7, q3, {2.983050737973751, -90.325884432736274},
                  {30.902915765072844, -6.396035046006810},
                  {-6.822717212205538, -1.551498996606969});
}

TEST(FormfactorBasic, Bipyramid4)
{
    double length = 10.;
    double height = 4;
    double height_ratio = .7;
    double alpha = 0.8;

    double tga = std::tan(alpha);
    double H2divLtga = height * 2. / length / tga;
    double ratioH2divLtga = height_ratio * height * 2. / length / tga;
    double volume = 1. / 6. * tga * length * length * length
                    * (2. - (1. - H2divLtga) * (1. - H2divLtga) * (1. - H2divLtga)
                       - (1. - ratioH2divLtga) * (1. - ratioH2divLtga) * (1. - ratioH2divLtga));

    Bipyramid4 particle(length, height, height_ratio, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    const double H = height * (1 + height_ratio);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), H * 1e-15);
    EXPECT_NEAR(H, particle.spanZ(&rot).hig(), H * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q3, q7, {71.667831484482932, -230.808201823553276},
                  {9.207591216136489, -11.771946031248554},
                  {153.865204621358430, -35.153360066115340});
}

TEST(FormfactorBasic, Cylinder)
{
    EXPECT_THROW(Cylinder(-3, -5), std::runtime_error);
    // TODO EXPECT_THROW(Cylinder(0, 5), std::runtime_error);

    double radius = 3.;
    double height = 5.;
    double volume = pi * radius * radius * height;

    Cylinder particle(radius, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());

    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    RotationX rectX(90 * Units::deg);
    EXPECT_NEAR(-radius, particle.spanZ(&rectX).low(), 1e-13);
    EXPECT_NEAR(+radius, particle.spanZ(&rectX).hig(), 1e-13);
    RotationY rectY(90 * Units::deg);
    EXPECT_NEAR(-radius, particle.spanZ(&rectY).low(), 1e-13);
    EXPECT_NEAR(+radius, particle.spanZ(&rectY).hig(), 1e-13);

    RotationY invY(180 * Units::deg);
    EXPECT_NEAR(-height, particle.spanZ(&invY).low(), 1e-13);
    EXPECT_NEAR(0, particle.spanZ(&invY).hig(), 1e-13);

    for (double gamma : {1.123, -2.34, 7.5, -9.})
        // 7.5deg is worst case for 24-vertex circle
        EXPECT_NEAR(-radius,
                    particle.spanZ(new RotationEuler(0, 90 * Units::deg, gamma * Units::deg)).low(),
                    3e-2); // TODO decrease epsilon after replacement of vertex-based approximation

    test_ff(&particle);
    test_ff_print(&particle, q8, q5, q1, {18.469449230629245, 5.713270157358191},
                  {17.176432808867649, -22.792734170623486},
                  {-4.645065790555677, -29.883784273353299});
}

TEST(FormfactorBasic, HorizontalCylinder)
{
    EXPECT_THROW(HorizontalCylinder(-3, -5), std::runtime_error);

    double R = 3.; // radius
    double L = 5.; // length

    double b = -0;  // slice_bottom
    double t = 2.3; // slice_top

    double full_volume = pi * R * R * L;
    double trunc_volume =
        L
        * (t * sqrt(R * R - t * t) - b * sqrt(R * R - b * b) + R * R * (asin(t / R) - asin(b / R)));

    RotationZ rot(.42);

    // full

    HorizontalCylinder particle_full(R, L);
    EXPECT_DOUBLE_EQ(full_volume, particle_full.volume());

    EXPECT_NEAR(0., particle_full.spanZ(&rot).low(), 2 * R * 1e-15);
    EXPECT_NEAR(2 * R, particle_full.spanZ(&rot).hig(), 2 * R * 1e-15);

    // truncated

    HorizontalCylinder particle_top(R, L, 0., R);
    EXPECT_NEAR(full_volume / 2, particle_top.volume(), 1e-7);

    HorizontalCylinder particle_bottom(R, L, -R, 0.);
    EXPECT_NEAR(full_volume / 2, particle_bottom.volume(), 1e-7);

    HorizontalCylinder particle_trunc(R, L, b, t);
    EXPECT_NEAR(trunc_volume, particle_trunc.volume(), 1e-7);

    EXPECT_NEAR(0., particle_trunc.spanZ(&rot).low(), 1e-7);
    EXPECT_NEAR(t - b, particle_trunc.spanZ(&rot).hig(), 1e-7);

    RotationX rot1((pi / 2));
    EXPECT_NEAR(-R, particle_trunc.spanZ(&rot1).low(), 1e-7);
    EXPECT_NEAR(R, particle_trunc.spanZ(&rot1).hig(), 1e-7);

    RotationX rot2((pi / 4));
    EXPECT_NEAR(-R / sqrt(2.), particle_trunc.spanZ(&rot2).low(), 1e-7);
    EXPECT_NEAR(R, particle_trunc.spanZ(&rot2).hig(), 1e-7);


    test_ff(&particle_trunc);
    test_ff_print(&particle_trunc, q5, q3, q8, {14.015342836573883, -6.888112822293706},
                  {14.242986207014122, 8.418849429673980}, {9.515192299073808, 1.605174270540351});
}

TEST(FormfactorBasic, Dodecahedron)
{
    EXPECT_THROW(Dodecahedron(-1), std::runtime_error);
    EXPECT_THROW(Dodecahedron(0), std::runtime_error);

    double edge = 3.;
    double volume = (15 + 7 * sqrt(5)) / 4 * pow(edge, 3);

    Dodecahedron particle(edge);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), 1e-4);
    EXPECT_NEAR(2 * 1.11352 * edge, particle.spanZ(&rot).hig(), 1e-4);
    // height=2*inradius from web ressource

    test_ff(&particle);
    test_ff_print(&particle, q1, q9, q2, {-10.878977464020377, -12.054560796788392},
                  {17.159983916914797, -101.323984000413375},
                  {3.473237908867695, -101.684458467045090});
}

TEST(FormfactorBasic, EllipsoidalCylinder)
{
    double radiusx = 3.;
    double radiusy = 5.;
    double height = 4;
    double volume = pi * radiusx * radiusy * height;

    EllipsoidalCylinder particle(radiusx, radiusy, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q3, q5, q7, {-10.885311020650205, -11.665937196046933},
                  {-10.799240867203938, 9.860673803151998},
                  {80.174590178312371, -9.667398944888935});
}

TEST(FormfactorBasic, Sphere)
{
    EXPECT_THROW(Sphere(-1), std::runtime_error);
    // TODO EXPECT_THROW(Sphere(0), std::runtime_error);

    double radius = 5.;
    double volume = 4. / 3. * pi * radius * radius * radius;

    Sphere particle(radius);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), 2 * radius * 1e-15);
    EXPECT_NEAR(2 * radius, particle.spanZ(&rot).hig(), 2 * radius * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q7, q2, q9, {244.485507285703591, -75.628229907665755},
                  {-78.931070547231570, -88.340730913140689},
                  {-61.244894114710341, -104.719370472590086});
}

TEST(FormfactorBasic, Spheroid)
{
    EXPECT_THROW(Spheroid(-3, -5), std::runtime_error);
    // TODO EXPECT_THROW(Spheroid(0, 5), std::runtime_error);
    // TODO EXPECT_THROW(Spheroid(3, 0), std::runtime_error);
    double radius = 3.;
    double height = 5.;
    double volume = 2. / 3. * pi * radius * radius * height;

    Spheroid particle(radius, height / 2);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q1, q4, q6, {-4.602530455511736, -29.610135452483558},
                  {32.081363863555183, -59.315827480130920},
                  {29.477269660229449, 23.594341882593568});
}

TEST(FormfactorBasic, HemiEllipsoid)
{
    EXPECT_THROW(HemiEllipsoid(-6, -7, 5), std::runtime_error);

    double radiusx = 6.;
    double radiusy = 7.;
    double height = 5.;

    double volume = (2 * pi) * radiusx * radiusy * height / 3.;

    HemiEllipsoid particle(radiusx, radiusy, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());

    test_ff(&particle);
}

TEST(FormfactorBasic, Icosahedron)
{
    EXPECT_THROW(Icosahedron(-1), std::runtime_error);
    EXPECT_THROW(Icosahedron(0), std::runtime_error);

    double edge = 7.;
    double volume = 5 * (3 + sqrt(5)) / 12 * pow(edge, 3);

    Icosahedron particle(edge);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), 1e-4);
    EXPECT_NEAR(2 * 0.755761 * edge, particle.spanZ(&rot).hig(), 1e-4);
    // height=2*inradius from web ressource

    test_ff(&particle);
    test_ff_print(&particle, q1, q8, q9, {59.253013256524476, -33.045312405031297},
                  {-65.498197697297627, -48.242910954959228},
                  {-56.385591749217589, -73.996078291078817});
}

TEST(FormfactorBasic, PlatonicOctahedron)
{
    EXPECT_THROW(PlatonicOctahedron(-1), std::runtime_error);
    EXPECT_THROW(PlatonicOctahedron(0), std::runtime_error);

    double edge = 7.;
    double volume = sqrt(2.) / 3 * pow(edge, 3);

    PlatonicOctahedron particle(edge);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), 1e-4);
    EXPECT_NEAR(edge * sqrt(2.), particle.spanZ(&rot).hig(), 1e-4);

    test_ff(&particle);
    test_ff_print(&particle, q3, q2, q7, {-14.326152086303487, 29.017548785469330},
                  {-55.305324166329513, -64.855380598892353},
                  {111.736751011333311, -34.195429173985289});
}

TEST(FormfactorBasic, PlatonicTetrahedron)
{
    EXPECT_THROW(PlatonicTetrahedron(-1), std::runtime_error);
    EXPECT_THROW(PlatonicTetrahedron(0), std::runtime_error);

    double edge = 7.;
    double volume = sqrt(2.) / 12 * pow(edge, 3);

    PlatonicTetrahedron particle(edge);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), 1e-4);
    EXPECT_NEAR(edge * sqrt(2. / 3.), particle.spanZ(&rot).hig(), 1e-4);

    test_ff(&particle);
    test_ff_print(&particle, q1, q0, q4, {7.338200374299324, -18.083987657244137},
                  {32.356617373497095, -16.489570069327176},
                  {25.457763901100478, -17.975973097181353});
}

TEST(FormfactorBasic, Prism3)
{
    EXPECT_THROW(Prism3(-4, -6), std::runtime_error);
    // TODO EXPECT_THROW(Prism3(0, 6), std::runtime_error);
    // TODO EXPECT_THROW(Prism3(4, 0), std::runtime_error);

    double height = 4.;
    double base_edge = 6.;
    double volume = sqrt(3.) / 4. * height * base_edge * base_edge;

    Prism3 particle(base_edge, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q5, q3, {43.549331118825300, -33.846003647829569},
                  {14.780314930820900, -21.332048794901649},
                  {22.325921825843015, 18.106569746121767});
}

TEST(FormfactorBasic, Prism6)
{
    EXPECT_THROW(Prism6(-4, -3), std::runtime_error);
    // TODO EXPECT_THROW(Prism6(0, 3), std::runtime_error);
    // TODO EXPECT_THROW(Prism6(4, 0), std::runtime_error);

    double height = 4.;
    double base_edge = 3.;
    double volume = 3. * sqrt(3.) / 2. * height * base_edge * base_edge;

    Prism6 particle(base_edge, height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q6, q4, {64.541920636979285, -50.091301690725452},
                  {30.516218708711495, 18.292325480907053},
                  {43.176189009775861, -50.151555255241675});
}

TEST(FormfactorBasic, Pyramid4)
{
    EXPECT_THROW(Pyramid4(-10, -4, .8), std::runtime_error);
    EXPECT_THROW(Pyramid4(0, 4, .8), std::runtime_error);
    // TODO EXPECT_THROW(Pyramid4(10, 0, .8), std::runtime_error);
    EXPECT_THROW(Pyramid4(10, 4, 0), std::runtime_error);

    double base_edge = 10.;
    double height = 4.;
    double alpha = 0.8;
    double tga = std::tan(alpha);
    double H2divLtga = height * 2. / base_edge / tga;
    double volume = 1. / 6. * tga * base_edge * base_edge * base_edge
                    * (1. - (1. - H2divLtga) * (1. - H2divLtga) * (1. - H2divLtga));

    Pyramid4 particle(base_edge, height, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q6, q8, q0, {-3.441841162284046, 8.991520427160811},
                  {4.963277061028782, 3.351425235955785},
                  {123.934890528173085, -57.248659449949685});
}

TEST(FormfactorBasic, SawtoothRippleBox)
{
    EXPECT_THROW(SawtoothRippleBox(-20, 4, -100, .3), std::runtime_error);
    // TODO EXPECT_THROW(SawtoothRippleBox(0, 4, 100, .3), std::runtime_error);
    // TODO EXPECT_THROW(SawtoothRippleBox(20, 0, 100, .3), std::runtime_error);
    // TODO EXPECT_THROW(SawtoothRippleBox(20, 4, 0, .3), std::runtime_error);
    // TODO EXPECT_THROW(SawtoothRippleBox(20, 4, 0, -.3), std::runtime_error);

    double width = 20.;
    double height = 4.;
    double length = 100.0;
    double d = 0.3; // asymmetry
    double volume = 0.5 * height * width * length;

    SawtoothRippleBox particle(length, width, height, d);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    // test_ff( &particle ); WAITING: restore once radius returns the umkreis radius
    test_ff_print(&particle, q9, q2, q6, {34.804843364383807, -22.650035449823580},
                  {-15.198443862633562, -71.676907701099211},
                  {-7.303609141120358, -5.907207608040455});
}

TEST(FormfactorBasic, TruncatedCube)
{
    EXPECT_THROW(TruncatedCube(0, 0), std::runtime_error);
    EXPECT_THROW(TruncatedCube(10, -1), std::runtime_error);
    EXPECT_THROW(TruncatedCube(10, 6), std::runtime_error);

    double length = 15.;
    double t = 6.; // side length of removed trirectangular tetrahedron at each vertex
    double volume = length * length * length - 4. / 3. * t * t * t;

    TruncatedCube particle(length, t);
    EXPECT_DOUBLE_EQ(particle.volume(), volume);
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), length * 1e-15);
    EXPECT_NEAR(length, particle.spanZ(&rot).hig(), length * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q1, q7, {-354.587043898980710, -278.961804308277806},
                  {3.203168895902472, 6.422974142279478},
                  {-130.107009837376836, 62.848850174170870});
}

TEST(FormfactorBasic, SphericalSegment)
{
    EXPECT_THROW(SphericalSegment(-5, -3, 1.7), std::runtime_error);
    // TODO EXPECT_THROW(SphericalSegment(5, 3, 0), std::runtime_error);
    EXPECT_THROW(SphericalSegment(0, 3, 1.7), std::runtime_error);
    EXPECT_THROW(SphericalSegment(5, 8.4, 1.7), std::runtime_error);

    double radius = 5.;
    double height = 3.;
    double HdivR = height / radius;
    double volume = pi / 3. * radius * radius * radius
                    * (3. * HdivR - 1. - (HdivR - 1.) * (HdivR - 1.) * (HdivR - 1.));

    SphericalSegment particle(radius, 0, 2 * radius - height);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q8, q5, q3, {-2.581178686741104, 1.096481712394733},
                  {1.598578167217787, -5.876368926542401}, {7.334079342958564, 9.540086846266590});
}

TEST(FormfactorBasic, SpheroidalSegment)
{
    EXPECT_THROW(SpheroidalSegment(-5, -3, 1.7, 0), std::runtime_error);
    // TODO EXPECT_THROW(SpheroidalSegment(5, 3, 0, 0), std::runtime_error);
    EXPECT_THROW(SpheroidalSegment(0, 3, 1.7, 0), std::runtime_error);
    // TODO EXPECT_THROW(SpheroidalSegment(5, 0, 1.7, 0), std::runtime_error);

    double radius = 3.;
    double height = 5.;
    double flattening = 1.5;
    double volume =
        pi * radius * height * height / flattening * (1. - height / (3. * flattening * radius));

    double radius_xy = radius;
    double radius_z = radius_xy * flattening;
    double rm_bottom = 2 * radius_z - height;

    SpheroidalSegment particle(radius_xy, radius_z, 0, rm_bottom);
    EXPECT_DOUBLE_EQ(volume, particle.volume());
    RotationZ rot(.42);
    EXPECT_NEAR(0., particle.spanZ(&rot).low(), height * 1e-15);
    EXPECT_NEAR(height, particle.spanZ(&rot).hig(), height * 1e-15);

    test_ff(&particle);
    test_ff_print(&particle, q0, q9, q5, {69.033474560055168, -50.648848240714358},
                  {43.816907833926741, -48.128263322981383},
                  {18.861884199850003, -22.594798796005833});
}

TEST(FormfactorBasic, Pyramid3)
{
    EXPECT_THROW(Pyramid3(-16, -4, .8), std::runtime_error);
    EXPECT_THROW(Pyramid3(0, 4, .8), std::runtime_error);
    // TODO EXPECT_THROW(Pyramid3(16, 0, .8), std::runtime_error);
    EXPECT_THROW(Pyramid3(16, 4, 0), std::runtime_error);

    double base_edge = 16.;
    double height = 4.;
    double alpha = 0.8;
    double tga = std::tan(alpha);
    double sqrt3H2divLtga = std::sqrt(3.) * 2. * height / base_edge / tga;
    double volume = tga / 24. * base_edge * base_edge * base_edge
                    * (1. - (1. - sqrt3H2divLtga) * (1. - sqrt3H2divLtga) * (1. - sqrt3H2divLtga));

    Pyramid3 particle(base_edge, height, alpha);
    EXPECT_DOUBLE_EQ(volume, particle.volume());

    test_ff(&particle);
    test_ff_print(&particle, q6, q8, q3, {12.744695481874178, -7.093565343710422},
                  {-5.320188956629416, -14.254410426450482},
                  {9.442516587003057, -7.290344882322833});
}