IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 528d16d2 authored by ralfulrich's avatar ralfulrich
Browse files

work on unit tests

parent 41ab23c7
No related branches found
No related tags found
1 merge request!369Resolve "Trajectories crossing observation plane"
Pipeline #4968 passed
...@@ -36,7 +36,7 @@ namespace corsika { ...@@ -36,7 +36,7 @@ namespace corsika {
std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon); std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
if (!x3.size()) { if (!x3.size()) {
return {}; // no solution, numeric problem return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
} }
long double y = x3[0]; // there is always at least one solution long double y = x3[0]; // there is always at least one solution
// The essence - choosing Y with maximal absolute value. // The essence - choosing Y with maximal absolute value.
......
...@@ -31,6 +31,7 @@ namespace corsika { ...@@ -31,6 +31,7 @@ namespace corsika {
// tracking: Note, this is NOT a general solution and should be clearly revised with // tracking: Note, this is NOT a general solution and should be clearly revised with
// a more robust tracking. #ifdef DEBUG // a more robust tracking. #ifdef DEBUG
if (deleteOnHit_) { if (deleteOnHit_) {
// since this is basically a bug, it cannot be tested LCOV_EXCL_START
LengthType const check = LengthType const check =
(particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal()); (particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal());
if (check < 0_m) { if (check < 0_m) {
...@@ -38,6 +39,7 @@ namespace corsika { ...@@ -38,6 +39,7 @@ namespace corsika {
CORSIKA_LOG_WARN("Temporary fix: write and remove particle."); CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
} else } else
return ProcessReturn::Ok; return ProcessReturn::Ok;
// LCOV_EXCL_STOP
} else } else
// #endif // #endif
return ProcessReturn::Ok; return ProcessReturn::Ok;
......
...@@ -32,7 +32,7 @@ namespace corsika { ...@@ -32,7 +32,7 @@ namespace corsika {
* $ * $
* \vec{x}(t_{i+0.5}) = \vec{x}(t_{i}) + \vec{v(t_{i})} * \Delta t / 2 \\ * \vec{x}(t_{i+0.5}) = \vec{x}(t_{i}) + \vec{v(t_{i})} * \Delta t / 2 \\
* \vec{v}(t_{i+1}) = \vec{v}(t_{i}) + \vec{v}(t_{i})\cross\vec{B}(x_{i}, t_{i}) * * \vec{v}(t_{i+1}) = \vec{v}(t_{i}) + \vec{v}(t_{i})\cross\vec{B}(x_{i}, t_{i}) *
*\Delta t \\ * \Delta t \\
* \vec{x}(t_{i+1}) = \vec{x}(t_{i+0.5}) + \vec{v}(t_{i+1}) * \Delta t /2 \\ * \vec{x}(t_{i+1}) = \vec{x}(t_{i+0.5}) + \vec{v}(t_{i+1}) * \Delta t /2 \\
* $ * $
* *
......
...@@ -110,7 +110,6 @@ private: ...@@ -110,7 +110,6 @@ private:
Plane plane_; Plane plane_;
}; };
template <typename T> template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
......
...@@ -305,10 +305,12 @@ TEST_CASE("Solver") { ...@@ -305,10 +305,12 @@ TEST_CASE("Solver") {
pol4(z4, a, b, c, d, e)); pol4(z4, a, b, c, d, e));
vector<double> s1 = andre::solve_quartic_real(a, b, c, d, e); vector<double> s1 = andre::solve_quartic_real(a, b, c, d, e);
// vector<double> s1 = solve_quartic_real(a, b, c, d, e, epsilon);
remove_duplicates(s1, epsilon_check * 10); remove_duplicates(s1, epsilon_check * 10);
vector<double> s2 = solve_quartic_real(a, b, c, d, e);
remove_duplicates(s2, epsilon_check * 5);
CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", ")); CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
CORSIKA_LOG_INFO("N={}, s2=[{}]", s2.size(), fmt::join(s2, ", "));
CHECK(s1.size() == idegree + 1); CHECK(s1.size() == idegree + 1);
for (double value : s1) { for (double value : s1) {
...@@ -326,8 +328,27 @@ TEST_CASE("Solver") { ...@@ -326,8 +328,27 @@ TEST_CASE("Solver") {
(value == Approx(z4).epsilon(epsilon_check)))); (value == Approx(z4).epsilon(epsilon_check))));
} }
} }
// this is a bit less precise
CHECK(s2.size() == idegree + 1);
for (double value : s2) {
CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} z4={} eps_check={}", value, z1,
z2, z3, z4, epsilon_check);
if (std::abs(value) < epsilon_check) {
CHECK(((value == Approx(z1).margin(epsilon_check * 5)) ||
(value == Approx(z2).margin(epsilon_check * 5)) ||
(value == Approx(z3).margin(epsilon_check * 5)) ||
(value == Approx(z4).margin(epsilon_check * 5))));
} else {
CHECK(((value == Approx(z1).epsilon(epsilon_check * 5)) ||
(value == Approx(z2).epsilon(epsilon_check * 5)) ||
(value == Approx(z3).epsilon(epsilon_check * 5)) ||
(value == Approx(z4).epsilon(epsilon_check * 5))));
}
}
} }
} }
} }
} // quartic } // quartic
} }
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