test_bulk_sphere.cc 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129
  1. #include <complex>
  2. #include "../src/mesomie.hpp"
  3. #include "../src/nmie-basic.hpp"
  4. #include "gtest/gtest.h"
  5. // TODO fails for MP with 100 digits. And 16 digits, which should be equal to
  6. // double precision.
  7. #ifndef MULTI_PRECISION
  8. // TEST(BulkSphere, DISABLED_ArgPi) {
  9. //******************************************************************************
  10. TEST(BulkSphere, ArgPi) {
  11. std::vector<double> WLs{50, 80, 100, 200, 400}; // nm
  12. double host_index = 2.;
  13. double core_radius = 100.; // nm
  14. double delta = 1e-5;
  15. nmie::MultiLayerMie<nmie::FloatType> nmie;
  16. nmie.SetLayersIndex({std::complex<double>(4, 0)});
  17. for (auto WL : WLs) {
  18. nmie.SetLayersSize(
  19. {2 * nmie::PI_ * host_index * core_radius / (WL + delta)});
  20. nmie.RunMieCalculation();
  21. double Qabs_p = std::abs(static_cast<double>(nmie.GetQabs()));
  22. nmie.SetLayersSize(
  23. {2 * nmie::PI_ * host_index * core_radius / (WL - delta)});
  24. nmie.RunMieCalculation();
  25. double Qabs_m = std::abs(static_cast<double>(nmie.GetQabs()));
  26. nmie.SetLayersSize({2 * nmie::PI_ * host_index * core_radius / (WL)});
  27. nmie.RunMieCalculation();
  28. double Qabs = std::abs(static_cast<double>(nmie.GetQabs()));
  29. EXPECT_GT(Qabs_p + Qabs_m, Qabs);
  30. }
  31. }
  32. #endif
  33. //******************************************************************************
  34. // A list of tests for a bulk sphere from
  35. // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
  36. // table 1: sphere size and refractive index
  37. // followed by resulting extinction and scattering efficiencies
  38. std::vector<std::tuple<double, std::complex<double>, double, double, char> >
  39. parameters_and_results{
  40. // x, {Re(m), Im(m)}, Qext, Qsca, test_name
  41. // {0.099, {0.75, 0}, 7.417859e-06, 7.417859e-06, 'a'},
  42. // {0.099, {1.75, 0}, 7.417859e-06, 7.417859e-06, 'z'},
  43. // {0.101, {0.75, 0}, 8.033538e-06, 8.033538e-06, 'b'},
  44. // {10, {0.75, 0}, 2.232265, 2.232265, 'c'},
  45. // {100, {1.33, 1e-5}, 2.101321, 2.096594, 'e'},
  46. // {0.055, {1.5, 1}, 0.10149104, 1.131687e-05, 'g'},
  47. // {0.056, {1.5, 1}, 0.1033467, 1.216311e-05, 'h'},
  48. // {100, {1.5, 1}, 2.097502, 1.283697, 'i'},
  49. // {1, {10, 10}, 2.532993, 2.049405, 'k'},
  50. // {1000, {0.75, 0}, 1.997908, 1.997908, 'd'},
  51. {100, {10, 10}, 2.071124, 1.836785, 'l'},
  52. // {10000, {1.33, 1e-5}, 2.004089, 1.723857, 'f'},
  53. // {10000, {1.5, 1}, 2.004368, 1.236574, 'j'},
  54. // {10000, {10, 10}, 2.005914, 1.795393, 'm'},
  55. };
  56. //******************************************************************************
  57. TEST(BulkSphere, DISABLED_HandlesInput) {
  58. // TEST(BulkSphere, HandlesInput) {
  59. nmie::MultiLayerMie<nmie::FloatType> nmie;
  60. for (const auto& data : parameters_and_results) {
  61. auto x = std::get<0>(data);
  62. auto m = std::get<1>(data);
  63. // auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
  64. nmie.SetLayersSize({x});
  65. nmie.SetLayersIndex({m});
  66. // nmie.SetMaxTerms(Nstop);
  67. nmie.RunMieCalculation();
  68. double Qext = static_cast<double>(nmie.GetQext());
  69. double Qsca = static_cast<double>(nmie.GetQsca());
  70. double Qext_Du = std::get<2>(data);
  71. double Qsca_Du = std::get<3>(data);
  72. EXPECT_FLOAT_EQ(Qext_Du, Qext)
  73. << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
  74. << "\nnmax_ = " << nmie.GetMaxTerms();
  75. EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
  76. << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
  77. }
  78. }
  79. //******************************************************************************
  80. TEST(BulkSphere, MesoMie) {
  81. nmie::MultiLayerMie<nmie::FloatType> nmie;
  82. nmie::MesoMie<nmie::FloatType> mesomie;
  83. for (const auto& data : parameters_and_results) {
  84. auto x = std::get<0>(data);
  85. auto m = std::get<1>(data);
  86. // auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
  87. nmie.SetLayersSize({x});
  88. nmie.SetLayersIndex({m});
  89. // nmie.SetMaxTerms(Nstop);
  90. nmie.calcScattCoeffs();
  91. auto an_nmie = nmie.GetAn();
  92. int nmax = an_nmie.size();
  93. mesomie.calc_ab(nmax + 2, // nmax
  94. x, // R
  95. {x, 0}, // xd
  96. x * m, // xm
  97. 1.0 / (x * m), // eps_d
  98. m / x, // eps_m
  99. {0, 0}, // d_parallel
  100. {0, 0}); // d_perp
  101. auto an_meso = mesomie.GetAn();
  102. mesomie.calc_ab_classic(nmax + 2, x, m);
  103. auto an_meso_cl = mesomie.GetAnClassic();
  104. for (int n = 0; n < nmax && n < 200; n++) {
  105. std::cout << n << an_nmie[n] << ' ' << an_meso[n + 1] << ' '
  106. << an_meso_cl[n + 1] << std::endl;
  107. }
  108. // double Qext = static_cast<double>(nmie.GetQext());
  109. // double Qsca = static_cast<double>(nmie.GetQsca());
  110. // double Qext_Du = std::get<2>(data);
  111. // double Qsca_Du = std::get<3>(data);
  112. // EXPECT_FLOAT_EQ(Qext_Du, Qext)
  113. // << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
  114. // << "\nnmax_ = " << nmie.GetMaxTerms();
  115. // EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
  116. // << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
  117. }
  118. }
  119. int main(int argc, char** argv) {
  120. testing::InitGoogleTest(&argc, argv);
  121. return RUN_ALL_TESTS();
  122. }