test_Riccati_Bessel_logarithmic_derivative.cc 6.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156
  1. #include "gtest/gtest.h"
  2. #include "../src/nmie-impl.hpp"
  3. #include "test_spec_functions_data.hpp"
  4. // From W. Yang APPLIED OPTICS Vol. 42, No. 9, 20 March 2003
  5. // Dtest refractive index is m={1.05,1}, the size parameter is x = 80
  6. std::vector<int> Dtest_n({0,1,30,50,60,70,75,80,85,90,99,116,130});
  7. std::vector< std::complex<double>>
  8. Dtest_D1({
  9. //Orig
  10. // {0.11449e-15 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 },
  11. // {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935},
  12. // {0.13645,-0.10019e+01 },{0.18439,-0.10070e+01 },
  13. // {0.21070,-0.10107e+01 },{0.23845,-0.10154e+01 },
  14. // {0.26752,-0.10210e+01 },{0.29777,-0.10278e+01 },
  15. // {0.35481,-0.10426e+01 },{0.46923,-0.10806e+01 },
  16. // {0.17656,-0.13895e+01 }
  17. // mod (from Python mpmath)
  18. {0.0,-1.0}, {7.464603828e-5,-0.9999958865},
  19. {0.03476380918,-0.9986960672},{0.09529213152,-0.999347654},
  20. {0.1364513887,-1.001895883},{0.184388335,-1.006979164},
  21. {0.2107044267,-1.01072099},{0.2384524295,-1.015382914},
  22. {0.2675164524,-1.021040337},{0.2977711192,-1.027753418},
  23. {0.3548096904,-1.042622957},{0.4692294405,-1.080629479},
  24. {0.5673827836,-1.121108944},
  25. });
  26. std::vector< std::complex<double>>
  27. Dtest_D2({{0.64966e-69 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 },
  28. {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935},
  29. {0.13645,-0.10019e+01 },{0.17769,-0.10099e+01 },
  30. {0.41264e-01 ,-0.21076e+01 },{-0.20190,0.10435e+01 },
  31. {-0.26343,0.10223e+01 },{-0.29339,0.10291e+01 },
  32. {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 },
  33. {-0.56047,0.11206e+01 }});
  34. std::vector< std::complex<double>>
  35. Dtest_D3({{0.00000,0.10000e+01 },{-0.73809e-04 ,0.10000e+01 },
  36. {-0.34344e-01 ,0.99912},{-0.94022e-01 ,0.10004e+01 },
  37. {-0.13455,0.10032e+01 },{-0.18172,0.10084e+01 },
  38. {-0.20762,0.10122e+01 },{-0.23494,0.10169e+01 },
  39. {-0.26357,0.10225e+01 },{-0.29339,0.10291e+01 },
  40. {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 },
  41. {-0.56047,0.11206e+01 }});
  42. int LeRu_cutoff(std::complex<double> z) {
  43. auto x = std::abs(z);
  44. return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
  45. }
  46. TEST(D1test, mpmath_generated_input) {
  47. double min_abs_tol = 2e-11;
  48. for (const auto &data : D1_test_16digits) {
  49. auto z = std::get<0>(data);
  50. auto n = std::get<1>(data);
  51. auto D1_mp = std::get<2>(data);
  52. auto re_abs_tol = ( std::get<3>(data) > min_abs_tol && std::real(D1_mp) < min_abs_tol)
  53. ? std::get<3>(data) : min_abs_tol;
  54. auto im_abs_tol = ( std::get<4>(data) > min_abs_tol && std::imag(D1_mp) < min_abs_tol)
  55. ? std::get<4>(data) : min_abs_tol;
  56. // if re(D1) < 0.5 then round will give 0. To avoid zero tolerance add one.
  57. re_abs_tol *= std::abs(std::round(std::real(D1_mp))) + 1;
  58. im_abs_tol *= std::abs(std::round(std::imag(D1_mp))) + 1;
  59. auto Nstop = LeRu_cutoff(z)+1;
  60. std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop+35), r;
  61. int valid_digits = 6;
  62. int nstar = nmie::getNStar(Nstop, z, valid_digits);
  63. r.resize(nstar);
  64. nmie::evalBackwardR(z,r);
  65. nmie::convertRtoD1(z, r, Db);
  66. if (n > Db.size()) continue;
  67. EXPECT_NEAR(std::real(Db[n]), std::real(D1_mp), re_abs_tol)
  68. << "Db at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
  69. EXPECT_NEAR(std::imag(Db[n]), std::imag(D1_mp), im_abs_tol)
  70. << "Db at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
  71. nmie::evalDownwardD1(z, Dold);
  72. if (n > Dold.size()) continue;
  73. EXPECT_NEAR(std::real(Dold[n]), std::real(D1_mp), re_abs_tol)
  74. << "Dold at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
  75. EXPECT_NEAR(std::imag(Dold[n]), std::imag(D1_mp), im_abs_tol)
  76. << "Dold at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
  77. }
  78. }
  79. //TEST(D1test, DISABLED_WYang_data){
  80. TEST(D1test, WYang_data){
  81. double abs_tol = 1e-9;
  82. int test_loss_digits = std::round(15 - std::log10(1/abs_tol));
  83. int Nstop = 131;
  84. std::vector<std::complex<nmie::FloatType>> Df(Nstop), Db(Nstop),Dold(Nstop), r;
  85. std::complex<nmie::FloatType> z(1.05,1);
  86. z = z*80.0;
  87. // eval D1 directly from backward recurrence
  88. nmie::evalDownwardD1(z, Dold);
  89. // eval forward recurrence
  90. r.resize(Nstop+1);
  91. nmie::evalForwardR(z, r);
  92. nmie::convertRtoD1(z, r, Df);
  93. // eval backward recurrence
  94. int valid_digits = 6;
  95. int nstar = nmie::getNStar(Nstop, z, valid_digits);
  96. r.resize(nstar);
  97. nmie::evalBackwardR(z,r);
  98. nmie::convertRtoD1(z, r, Db);
  99. for (int i = 0; i < Dtest_n.size(); i++) {
  100. int n = Dtest_n[i];
  101. int forward_loss_digits = nmie::evalKapteynNumberOfLostSignificantDigits(n, z);
  102. forward_loss_digits += 3; // Kapteyn is too optimistic
  103. if (test_loss_digits > forward_loss_digits ) {
  104. EXPECT_NEAR(std::real(Df[n]), std::real(Dtest_D1[i]),
  105. abs_tol) << "f at n=" << n << " lost digits = " << forward_loss_digits;
  106. EXPECT_NEAR(std::imag(Df[n]), std::imag(Dtest_D1[i]),
  107. abs_tol) << "f at n=" << n << " lost digits = " << forward_loss_digits;
  108. }
  109. EXPECT_NEAR(std::real(Db[n]), std::real(Dtest_D1[i]),
  110. abs_tol) << "b at n=" << n;
  111. EXPECT_NEAR(std::imag(Db[n]), std::imag(Dtest_D1[i]),
  112. abs_tol) << "b at n=" << n;
  113. if (n < Dold.size()-15) {
  114. EXPECT_NEAR(std::real(Dold[n]), std::real(Dtest_D1[i]),
  115. abs_tol) << "old at n=" << n;
  116. EXPECT_NEAR(std::imag(Dold[n]), std::imag(Dtest_D1[i]),
  117. abs_tol) << "old at n=" << n;
  118. }
  119. }
  120. }
  121. TEST(KaptyenTest, HandlesInput) {
  122. // H.Du APPLIED OPTICS, Vol. 43, No. 9, 20 March 2004
  123. double l = nmie::evalKapteynNumberOfLostSignificantDigits(80, std::complex<double>(100,100));
  124. EXPECT_EQ(l, 7)<<"Should be equal";
  125. std::complex<double> z(10000,0);
  126. l = nmie::evalKapteynNumberOfLostSignificantDigits(5070, z);
  127. EXPECT_EQ(l, 0)<<"Should be equal";
  128. // find NStar such that l_nstar(z) - l_nmax(z) >= valid_digits
  129. int NStar = nmie::getNStar(5070, z,6);
  130. EXPECT_GE(NStar, 10130);
  131. // const double pi=3.14159265358979323846;
  132. // z = std::complex<double>(100,100);
  133. // l = nmie::evalKapteynNumberOfLostSignificantDigits(1, z);
  134. // EXPECT_EQ(l, 0)<<"Should be equal";
  135. }
  136. int main(int argc, char **argv) {
  137. testing::InitGoogleTest(&argc, argv);
  138. return RUN_ALL_TESTS();
  139. }