Line data Source code
1 : /* 2 : * Copyright (c) 2017 Sebastian Weber, Henri Menke. All rights reserved. 3 : * 4 : * This file is part of the pairinteraction library. 5 : * 6 : * The pairinteraction library is free software: you can redistribute it and/or modify 7 : * it under the terms of the GNU Lesser General Public License as published by 8 : * the Free Software Foundation, either version 3 of the License, or 9 : * (at your option) any later version. 10 : * 11 : * The pairinteraction library is distributed in the hope that it will be useful, 12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 : * GNU Lesser General Public License for more details. 15 : * 16 : * You should have received a copy of the GNU Lesser General Public License 17 : * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>. 18 : */ 19 : 20 : #include "QuantumDefect.hpp" 21 : #include "SQLite.hpp" 22 : #include "Wavefunction.hpp" 23 : 24 : #define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN 25 : #include <doctest/doctest.h> 26 : 27 : #include <iostream> 28 : 29 : template <int l> 30 : struct Fixture { 31 : QuantumDefect qd; 32 10 : Fixture() : qd("Rb", 79, l, 1.5){}; 33 : }; 34 : 35 5 : TEST_CASE_TEMPLATE("model_potentials", T, Fixture<1>, Fixture<2>) // NOLINT 36 : { 37 4 : T const fixture; 38 4 : auto const qd = fixture.qd; 39 : // There could be better coverage 40 2 : CHECK(std::isnan(model_potential::V(qd, 0))); 41 2 : CHECK(std::isnan(model_potential::g(qd, 0))); 42 2 : } 43 : 44 5 : TEST_CASE_TEMPLATE("numerovs_method", T, Fixture<1>, Fixture<2>) // NOLINT 45 : { 46 4 : T const fixture; 47 4 : auto const qd = fixture.qd; 48 4 : Numerov N(qd); 49 4 : auto const &xy = N.integrate(); 50 : 51 : // Check for correct number of integration steps 52 2 : CHECK(xy.rows() == 12087); 53 : 54 : // Check for correct upper bound and decay to zero 55 2 : CHECK(xy(xy.rows() - 1, 0) <= std::sqrt(2 * qd.n * (qd.n + 15))); 56 2 : CHECK(xy(xy.rows() - 1, 1) == doctest::Approx(0.0).epsilon(1e-6)); 57 2 : } 58 : 59 : #ifdef WITH_GSL 60 5 : TEST_CASE_TEMPLATE("coulomb_functions", T, Fixture<1>, Fixture<2>) // NOLINT 61 : { 62 4 : T const fixture; 63 4 : auto const qd = fixture.qd; 64 4 : Whittaker W(qd); 65 4 : auto const &xy = W.integrate(); 66 : 67 : // Check for correct number of integration steps 68 2 : CHECK(xy.rows() == 12087); 69 : 70 : // Check for correct upper bound and decay to zero 71 2 : CHECK(xy(xy.rows() - 1, 0) <= 2 * qd.n * (qd.n + 15)); 72 2 : CHECK(xy(xy.rows() - 1, 1) == doctest::Approx(0.0).epsilon(1e-6)); 73 2 : } 74 : 75 5 : TEST_CASE_TEMPLATE("method_comparison", T, Fixture<1>, Fixture<2>) // NOLINT 76 : { 77 4 : T const fixture; 78 4 : auto const qd = fixture.qd; 79 4 : Numerov N(qd); 80 4 : Whittaker W(qd); 81 4 : auto const &nxy = N.integrate(); 82 4 : auto const &wxy = W.integrate(); 83 : 84 : // Check whether both have the same number of points 85 2 : CHECK(nxy.rows() == wxy.rows()); 86 2 : size_t n = nxy.rows(); 87 : 88 : // Compare pointwise 89 24176 : for (size_t i = 0; i < n; ++i) { 90 24174 : CHECK(std::sqrt(nxy(i, 0)) * nxy(i, 1) - wxy(i, 1) == doctest::Approx(0.0).epsilon(1e-2)); 91 : } 92 2 : } 93 5 : TEST_CASE_TEMPLATE("integration", T, Fixture<1>, Fixture<2>) // NOLINT 94 : { 95 4 : T const fixture; 96 4 : auto const qd = fixture.qd; 97 2 : double mu_n = IntegrateRadialElement<Numerov>(qd, 1, qd); 98 2 : double mu_w = IntegrateRadialElement<Whittaker>(qd, 1, qd); 99 2 : CHECK(mu_n == doctest::Approx(mu_w).scale(1e-3)); // corresponds to 0.1% deviation 100 2 : } 101 : #endif