#include "gtest/gtest.h" #include "../src/nmie-impl.hpp" #include "test_spec_functions_data.hpp" // From W. Yang APPLIED OPTICS Vol. 42, No. 9, 20 March 2003 // Dtest refractive index is m={1.05,1}, the size parameter is x = 80 std::vector Dtest_n({0,1,30,50,60,70,75,80,85,90,99,116,130}); std::vector< std::complex> Dtest_D1({ //Orig // {0.11449e-15 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 }, // {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935}, // {0.13645,-0.10019e+01 },{0.18439,-0.10070e+01 }, // {0.21070,-0.10107e+01 },{0.23845,-0.10154e+01 }, // {0.26752,-0.10210e+01 },{0.29777,-0.10278e+01 }, // {0.35481,-0.10426e+01 },{0.46923,-0.10806e+01 }, // {0.17656,-0.13895e+01 } // mod (from Python mpmath) {0.0,-1.0}, {7.464603828e-5,-0.9999958865}, {0.03476380918,-0.9986960672},{0.09529213152,-0.999347654}, {0.1364513887,-1.001895883},{0.184388335,-1.006979164}, {0.2107044267,-1.01072099},{0.2384524295,-1.015382914}, {0.2675164524,-1.021040337},{0.2977711192,-1.027753418}, {0.3548096904,-1.042622957},{0.4692294405,-1.080629479}, {0.5673827836,-1.121108944}, }); std::vector< std::complex> Dtest_D2({{0.64966e-69 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 }, {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935}, {0.13645,-0.10019e+01 },{0.17769,-0.10099e+01 }, {0.41264e-01 ,-0.21076e+01 },{-0.20190,0.10435e+01 }, {-0.26343,0.10223e+01 },{-0.29339,0.10291e+01 }, {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 }, {-0.56047,0.11206e+01 }}); std::vector< std::complex> Dtest_D3({{0.00000,0.10000e+01 },{-0.73809e-04 ,0.10000e+01 }, {-0.34344e-01 ,0.99912},{-0.94022e-01 ,0.10004e+01 }, {-0.13455,0.10032e+01 },{-0.18172,0.10084e+01 }, {-0.20762,0.10122e+01 },{-0.23494,0.10169e+01 }, {-0.26357,0.10225e+01 },{-0.29339,0.10291e+01 }, {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 }, {-0.56047,0.11206e+01 }}); int LeRu_cutoff(std::complex z) { auto x = std::abs(z); return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1); } void parse_mpmath_data(const double min_abs_tol, const std::tuple< std::complex, int, std::complex, double, double > data, std::complex &z, int &n, std::complex &func_mp, double &re_abs_tol, double &im_abs_tol){ z = std::get<0>(data); n = std::get<1>(data); func_mp = std::get<2>(data); re_abs_tol = ( std::get<3>(data) > min_abs_tol && std::real(func_mp) < min_abs_tol) ? std::get<3>(data) : min_abs_tol; im_abs_tol = ( std::get<4>(data) > min_abs_tol && std::imag(func_mp) < min_abs_tol) ? std::get<4>(data) : min_abs_tol; // if re(func_mp) < 0.5 then round will give 0. To avoid zero tolerance add one. re_abs_tol *= std::abs(std::round(std::real(func_mp))) + 1; im_abs_tol *= std::abs(std::round(std::imag(func_mp))) + 1; } void parse2_mpmath_data(const nmie::FloatType min_abs_tol, const std::tuple< nmie::FloatType, std::complex, int, std::complex, nmie::FloatType, nmie::FloatType > data, nmie::FloatType &x, std::complex &m, int &n, std::complex &func_mp, nmie::FloatType &re_abs_tol, nmie::FloatType &im_abs_tol){ x = std::get<0>(data); m = std::get<1>(data); n = std::get<2>(data); func_mp = std::get<3>(data); re_abs_tol = ( std::get<4>(data) > min_abs_tol && std::real(func_mp) < min_abs_tol) ? std::get<4>(data) : min_abs_tol; im_abs_tol = ( std::get<5>(data) > min_abs_tol && std::imag(func_mp) < min_abs_tol) ? std::get<5>(data) : min_abs_tol; // if re(func_mp) < 0.5 then round will give 0. To avoid zero tolerance add one. re_abs_tol *= std::abs(std::round(std::real(func_mp))) + 1; im_abs_tol *= std::abs(std::round(std::imag(func_mp))) + 1; } template inline T pow2(const T value) {return value*value;} TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) { //TEST(zeta_psizeta_test, mpmath_generated_input) { double min_abs_tol = 2e-10; std::complex z, zeta_mp; int n; double re_abs_tol, im_abs_tol; for (const auto &data : zeta_test_16digits) { parse_mpmath_data(min_abs_tol, data, z, n, zeta_mp, re_abs_tol, im_abs_tol); auto Nstop = LeRu_cutoff(z)+10000; if (n > Nstop) continue; std::vector> D1dr(Nstop+135), D3(Nstop+135), PsiZeta(Nstop+135), Psi(Nstop); nmie::evalDownwardD1(z, D1dr); nmie::evalUpwardD3(z, D1dr, D3, PsiZeta); nmie::evalUpwardPsi(z, D1dr, Psi); auto a = std::real(PsiZeta[n]); auto b = std::imag(PsiZeta[n]); auto c = std::real(Psi[n]); auto d = std::imag(Psi[n]); auto c_one = std::complex(0, 1); auto zeta = (a*c + b*d)/(pow2(c) + pow2(d)) + c_one * ((b*c - a*d)/(pow2(c) + pow2(d))); // zeta = PsiZeta[n]/Psi[n]; if (std::isnan(std::real(zeta)) || std::isnan(std::imag(zeta))) continue; // std::vector> D1dr(Nstop+35), D3(Nstop+35), zeta(Nstop); // nmie::evalDownwardD1(z, D1dr); // nmie::evalUpwardD3(z, D1dr, D3); // nmie::evalUpwardZeta(z, D3, zeta); EXPECT_NEAR(std::real(zeta), std::real(zeta_mp), re_abs_tol) << "zeta at n=" << n << " Nstop="<< Nstop<<" z="< z(10000,0); l = nmie::evalKapteynNumberOfLostSignificantDigits(5070, z); EXPECT_EQ(l, 0)<<"Should be equal"; // find NStar such that l_nstar(z) - l_nmax(z) >= valid_digits int NStar = nmie::getNStar(5070, z,6); EXPECT_GE(NStar, 10130); // const double pi=3.14159265358979323846; // z = std::complex(100,100); // l = nmie::evalKapteynNumberOfLostSignificantDigits(1, z); // EXPECT_EQ(l, 0)<<"Should be equal"; } int main(int argc, char **argv) { testing::InitGoogleTest(&argc, argv); return RUN_ALL_TESTS(); }