1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
#ifdef ALGORITHM_DIAGNOSTIC
#include "Sample/HardParticle/HardParticles.h"
#include "Tests/GTestWrapper/google_test.h"
#include <cassert>
#include <complex>
#include <ff/Face.h> // ??
#include <iomanip>
#include <iostream>
#include <numbers>
#include <vector>
using std::numbers::pi;
const auto qlist = testing::Combine(
testing::Values(C3({1, 0, 0}), C3({0, 1, 0}), C3({0, 0, 1}), C3({1, 1, 0}), C3({1, 0, 1}),
C3({0, 1, 1}), C3({2, 3, 0}), C3({5, 0, 2}), C3({0, 5, 3}), C3({1, sqrt(2), 0}),
C3({sqrt(3), 0, 1}), C3({1, 1, 1})),
testing::Values(C3({1, 0, 0}), C3({0, 1, 0}), C3({0, 0, 1}), C3({1, 1, 0}), C3({1, 0, 1}),
C3({0, 1, 1}), C3({2, 3, 0}), C3({5, 0, 2}), C3({0, 5, 3}), C3({1, sqrt(2), 0}),
C3({sqrt(3), 0, 1}), C3({1, 1, 1})),
testing::Values(1e-19, 1e-17, 1e-15, 1e-12, 1e-9, 1e-6, 1e-3, .03, 1., 3., 10., 30., 100.),
testing::Values(-1., 1.), testing::Values(0, -4e-16, +8e-16, -5e-11, 3e-7, -2e-3, .01, .1));
complex_t deriv(const IFormfactor& ff, const C3& qf, const complex_t Ff,
const PolyhedralDiagnosis& /*df*/, const C3& qdir, const double qstep)
{
assert(Ff == ff.formfactor(qf));
complex_t Fi = ff.formfactor(qf + qstep * qdir);
PolyhedralDiagnosis di = polyhedralDiagnosis;
// assert(di==df);
return (Ff - Fi) / qstep;
}
//! Bisect between two q's to find possible discontinuities
void bisect(int& ifail, const C3& qdir, const IFormfactor& ff, const C3& q0, const C3& q1,
const complex_t F0, const complex_t F1, const PolyhedralDiagnosis& d0,
const PolyhedralDiagnosis& d1, const double qmindiff, const double Fmaxreldiff)
{
assert(d0 != d1);
if ((q0 - q1).mag() < qmindiff) {
// narrowed down to minimal step, now check for continuity
double aval = (std::abs(F0) + std::abs(F1)) / 2;
double step = std::abs(F0 - F1);
double relstep = step / aval;
if (relstep > Fmaxreldiff) {
std::cout << d0.message() << " -> " << d1.message() << ":\n";
std::cout << "relstep " << std::setprecision(8) << relstep << "=" << step << "/"
<< std::setprecision(16) << aval << "\n";
std::cout << " q[-] = " << q0 << "\n";
std::cout << " q[+] = " << q1 << "\n";
std::cout << " F[-] = " << F0 << "\n";
std::cout << " F[+] = " << F1 << "\n";
std::cout << " F'[-1k] =" << -deriv(ff, q0, F0, d0, -qdir, 1000 * qmindiff) << "\n";
std::cout << " F'[-300] =" << -deriv(ff, q0, F0, d0, -qdir, 300 * qmindiff) << "\n";
std::cout << " F'[-100] =" << -deriv(ff, q0, F0, d0, -qdir, 100 * qmindiff) << "\n";
std::cout << " F'[-30] =" << -deriv(ff, q0, F0, d0, -qdir, 30 * qmindiff) << "\n";
std::cout << " F'[-10] =" << -deriv(ff, q0, F0, d0, -qdir, 10 * qmindiff) << "\n";
std::cout << " F'[-3] =" << -deriv(ff, q0, F0, d0, -qdir, 3 * qmindiff) << "\n";
std::cout << " F'[-1] =" << -deriv(ff, q0, F0, d0, -qdir, 1 * qmindiff) << "\n";
std::cout << " F'[here]=" << (F1 - F0) / (q0 - q1).mag() << "\n";
std::cout << " F'[+1] =" << deriv(ff, q1, F1, d1, +qdir, 1 * qmindiff) << "\n";
std::cout << " F'[+3] =" << deriv(ff, q1, F1, d1, +qdir, 3 * qmindiff) << "\n";
std::cout << " F'[+10] =" << deriv(ff, q1, F1, d1, +qdir, 10 * qmindiff) << "\n";
std::cout << " F'[+30] =" << deriv(ff, q1, F1, d1, +qdir, 30 * qmindiff) << "\n";
std::cout << " F'[+100] =" << deriv(ff, q1, F1, d1, +qdir, 100 * qmindiff) << "\n";
std::cout << " F'[+300] =" << deriv(ff, q1, F1, d1, +qdir, 300 * qmindiff) << "\n";
std::cout << " F'[+1k] =" << deriv(ff, q1, F1, d1, +qdir, 1000 * qmindiff) << "\n";
std::cout << std::endl;
++ifail;
// maxrelstep = relstep;
return;
}
// std::cout<<"ok for "<<d0.message()<<"->"<<d1.message()<<" at q between "<<q0<<" and
// "<<q1<<std::endl;
return;
}
C3 q2 = (q0 + q1) / 2.;
complex_t F2 = ff.formfactor(q2);
PolyhedralDiagnosis d2 = polyhedralDiagnosis;
if (d2 != d0)
bisect(ifail, qdir, ff, q0, q2, F0, F2, d0, d2, qmindiff, Fmaxreldiff);
if (d2 != d1)
bisect(ifail, qdir, ff, q2, q1, F2, F1, d2, d1, qmindiff, Fmaxreldiff);
}
void run_bisection(int& ifail, IFormfactor& ff, const C3& q0, const C3& q1)
{
const double qdiffmin = std::max(q0.mag(), q1.mag()) / 4e11;
complex_t F0 = ff.formfactor(q0);
PolyhedralDiagnosis d0 = polyhedralDiagnosis;
complex_t F1 = ff.formfactor(q1);
PolyhedralDiagnosis d1 = polyhedralDiagnosis;
if (d0 == d1)
return;
bisect(ifail, q1 - q0, ff, q0, q1, F0, F1, d0, d1, qdiffmin, .6e-10);
}
void run_test(IFormfactor& ff)
{
::testing::internal::ParamGenerator<std::tuple<C3, C3, double, double, double>> gen = qlist;
int ifail = 0;
for (auto it : gen) {
const C3 q_dir0 = std::get<0>(it).unit_or_throw();
const C3 q_dir1 = std::get<1>(it).unit_or_throw();
const double qrealmag = std::get<2>(it);
const double qrel1 = std::get<3>(it);
const double qimagrel = std::get<4>(it);
const complex_t qmag(qrealmag, qrealmag * qimagrel);
const C3 q0 = q_dir0 * qmag;
const C3 q1 = q_dir1 * (qmag * qrel1);
run_bisection(ifail, ff, q0, q1);
}
EXPECT_EQ(ifail, 0);
}
class BisectFF : public testing::Test {};
TEST_F(BisectFF, Bisect1)
{
Pyramid6 ff(.23, 3.5, .999 * pi / 2);
// Pyramid3 ff(8.23, .27, .51);
run_test(ff);
}
#endif // ALGORITHM_DIAGNOSTIC
|
