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