IAP GITLAB
Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
corsika
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Issue analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Air Shower Physics
corsika
Commits
2c82a32f
Commit
2c82a32f
authored
3 years ago
by
ralfulrich
Committed by
Ralf Ulrich
3 years ago
Browse files
Options
Downloads
Patches
Plain Diff
coverage
parent
cc5279f0
No related branches found
Branches containing commit
No related tags found
Tags containing commit
1 merge request
!369
Resolve "Trajectories crossing observation plane"
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
corsika/detail/framework/utility/CubicSolver.inl
+4
-0
4 additions, 0 deletions
corsika/detail/framework/utility/CubicSolver.inl
tests/framework/testSolver.cpp
+94
-40
94 additions, 40 deletions
tests/framework/testSolver.cpp
with
98 additions
and
40 deletions
corsika/detail/framework/utility/CubicSolver.inl
+
4
−
0
View file @
2c82a32f
...
@@ -261,6 +261,9 @@ namespace corsika {
...
@@ -261,6 +261,9 @@ namespace corsika {
}
}
#endif
#endif
}
else
{
}
else
{
// this is only if the former solve_cubic_real_analytic would not result
// in any solution. We have no test case for this. This is excluded from tests:
// LCOV_EXCL_START
long
double
const
dist
=
std
::
fma
(
b
/
a
,
b
/
a
,
-
3
*
c
/
a
);
long
double
const
dist
=
std
::
fma
(
b
/
a
,
b
/
a
,
-
3
*
c
/
a
);
long
double
const
xinfl
=
-
b
/
(
a
*
3
);
long
double
const
xinfl
=
-
b
/
(
a
*
3
);
...
@@ -277,6 +280,7 @@ namespace corsika {
...
@@ -277,6 +280,7 @@ namespace corsika {
x1
=
xinfl
+
2
/
3
*
std
::
sqrt
(
dist
);
x1
=
xinfl
+
2
/
3
*
std
::
sqrt
(
dist
);
}
}
}
}
// LCOV_EXCL_STOP
}
}
long
double
f_x1
=
cubic_function
(
x1
,
a
,
b
,
c
,
d
);
long
double
f_x1
=
cubic_function
(
x1
,
a
,
b
,
c
,
d
);
...
...
This diff is collapsed.
Click to expand it.
tests/framework/testSolver.cpp
+
94
−
40
View file @
2c82a32f
...
@@ -98,35 +98,64 @@ TEST_CASE("Solver") {
...
@@ -98,35 +98,64 @@ TEST_CASE("Solver") {
long
double
b
=
-
z1
-
z2
;
long
double
b
=
-
z1
-
z2
;
long
double
c
=
z1
*
z2
;
long
double
c
=
z1
*
z2
;
std
::
vector
<
double
>
s1
=
solve_quadratic_real
(
a
,
b
,
c
);
{
remove_duplicates
(
s1
,
epsilon_check
);
std
::
vector
<
double
>
s1
=
solve_quadratic_real
(
a
,
b
,
c
);
remove_duplicates
(
s1
,
epsilon_check
);
CORSIKA_LOG_INFO
(
"quadratic: z1={}, z2={}, N={}, s1=[{}]"
,
z1
,
z2
,
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
CORSIKA_LOG_INFO
(
"quadratic: z1={}, z2={}, N={}, s1=[{}]"
,
z1
,
z2
,
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
CHECK
(
s1
.
size
()
==
idegree
+
1
);
for
(
auto
value
:
s1
)
{
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))));
}
else
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))));
}
}
}
CHECK
(
s1
.
size
()
==
idegree
+
1
);
// cubic with x^3 * 0 should result in the same
for
(
auto
value
:
s1
)
{
{
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
std
::
vector
<
double
>
s1
=
solve_cubic_real
(
0
,
a
,
b
,
c
);
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
remove_duplicates
(
s1
,
epsilon_check
);
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))));
}
else
{
CORSIKA_LOG_INFO
(
"quadratic/cubic: z1={}, z2={}, N={}, s1=[{}]"
,
z1
,
z2
,
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))));
CHECK
(
s1
.
size
()
==
idegree
+
1
);
for
(
auto
value
:
s1
)
{
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))));
}
else
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))));
}
}
}
}
}
/*
std::vector<complex<double>> s2 = solve_quadratic(a, b, c, epsilon);
// quartic with x^4 * 0 + x^3 * 0 must be the same
CHECK(s2.size() == idegree + 1);
{
for (auto value : s2) {
std
::
vector
<
double
>
s1
=
solve_quartic_real
(
0
,
0
,
a
,
b
,
c
);
if (std::abs(value) < epsilon) {
remove_duplicates
(
s1
,
epsilon_check
);
CHECK(((value == Approx(z1).margin(epsilon_check)) ||
(value == Approx(z2).margin(epsilon_check))));
CORSIKA_LOG_INFO
(
"quadratic/quartic: z1={}, z2={}, N={}, s1=[{}]"
,
z1
,
z2
,
} else {
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
(value == Approx(z2).epsilon(epsilon_check))));
CHECK
(
s1
.
size
()
==
idegree
+
1
);
for
(
auto
value
:
s1
)
{
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))));
}
else
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))));
}
}
}
}
}
}*/
}
}
}
}
}
}
...
@@ -167,23 +196,48 @@ TEST_CASE("Solver") {
...
@@ -167,23 +196,48 @@ TEST_CASE("Solver") {
"a={}, b={}, c={}, d={}"
,
"a={}, b={}, c={}, d={}"
,
z1
,
z2
,
z3
,
a
,
b
,
c
,
d
);
z1
,
z2
,
z3
,
a
,
b
,
c
,
d
);
vector
<
double
>
s1
=
solve_cubic_real
(
a
,
b
,
c
,
d
);
{
remove_duplicates
(
s1
,
epsilon_check
*
10
);
vector
<
double
>
s1
=
solve_cubic_real
(
a
,
b
,
c
,
d
);
remove_duplicates
(
s1
,
epsilon_check
*
10
);
CORSIKA_LOG_INFO
(
"N={}, s1=[{}]"
,
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
CORSIKA_LOG_INFO
(
"N={}, s1=[{}]"
,
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
CHECK
(
s1
.
size
()
==
idegree
+
1
);
for
(
double
value
:
s1
)
{
CORSIKA_LOG_INFO
(
"value={}, z1={} z2={} z3={} eps_check={}"
,
value
,
z1
,
z2
,
z3
,
epsilon_check
);
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z3
).
margin
(
epsilon_check
))));
}
else
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z3
).
epsilon
(
epsilon_check
))));
}
}
}
CHECK
(
s1
.
size
()
==
idegree
+
1
);
// quartic with x^4 *0 must be the same
for
(
double
value
:
s1
)
{
{
CORSIKA_LOG_INFO
(
"value={}, z1={} z2={} z3={} eps_check={}"
,
value
,
z1
,
z2
,
vector
<
double
>
s1
=
solve_quartic_real
(
0
,
a
,
b
,
c
,
d
);
z3
,
epsilon_check
);
remove_duplicates
(
s1
,
epsilon_check
*
10
);
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
CORSIKA_LOG_INFO
(
"(quartic) N={}, s1=[{}]"
,
s1
.
size
(),
fmt
::
join
(
s1
,
", "
));
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z3
).
margin
(
epsilon_check
))));
CHECK
(
s1
.
size
()
==
idegree
+
1
);
}
else
{
for
(
double
value
:
s1
)
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
CORSIKA_LOG_INFO
(
"value={}, z1={} z2={} z3={} eps_check={}"
,
value
,
z1
,
z2
,
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))
||
z3
,
epsilon_check
);
(
value
==
Approx
(
z3
).
epsilon
(
epsilon_check
))));
if
(
std
::
abs
(
value
)
<
epsilon_check
)
{
CHECK
(((
value
==
Approx
(
z1
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
margin
(
epsilon_check
))
||
(
value
==
Approx
(
z3
).
margin
(
epsilon_check
))));
}
else
{
CHECK
(((
value
==
Approx
(
z1
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z2
).
epsilon
(
epsilon_check
))
||
(
value
==
Approx
(
z3
).
epsilon
(
epsilon_check
))));
}
}
}
}
}
}
}
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment