Fix Nstar evaluation for z=0

src/nmie-basic.hpp

@@ -275,7 +275,7 @@ namespace nmie {
   // ********************************************************************** //
 
   template <typename FloatType>
-  unsigned int LeRu_cutoff(const std::complex<FloatType> zz) {
+  unsigned int LeRu_near_field_cutoff(const std::complex<FloatType> zz) {
     std::complex<double> z = ConvertComplex<double>(zz);
     auto x = std::abs(z);
     return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
@@ -297,8 +297,8 @@ namespace nmie {
     } else {
       nmax = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
     }
-    //Le Ru
-    auto Nstop = nmie::LeRu_cutoff(std::complex<FloatType>(xL,0))+1;
+    //Use Le Ru cutoff for near field, as a universal one.
+    auto Nstop = nmie::LeRu_near_field_cutoff(std::complex<FloatType>(xL, 0))+1;
     if (Nstop > nmax) nmax = Nstop;
     return nmax;
   }

src/nmie-nearfield.hpp

@@ -467,7 +467,8 @@ void MultiLayerMie<FloatType>::calcRadialOnlyDependantFunctions(const FloatType
   double delta_Rho = eval_delta<FloatType>(radius_points, from_Rho, to_Rho);
   for (unsigned int j=0; j < radius_points; j++) {
     auto Rho = static_cast<FloatType>(from_Rho + j*delta_Rho);
-    int near_field_nmax = calcNmax(to_Rho);
+//    if (Rho < 1e-5) Rho = 1e-5; // TODO do we need this?.
+    int near_field_nmax = calcNmax(Rho);
     // Skip if not enough terms in Mie series (i.e. required near field nmax > available terms )
     if (near_field_nmax > available_maximal_nmax_ && !isIgnoreAvailableNmax) continue;
     if (near_field_nmax > available_maximal_nmax_)  near_field_nmax = available_maximal_nmax_;

src/special-functions-impl.hpp

@@ -63,13 +63,15 @@ int evalKapteynNumberOfLostSignificantDigits(const int ni,
                                              const std::complex<FloatType> zz) {
   using std::abs;  using std::imag; using std::real; using std::log; using std::sqrt; using std::round;
   std::complex<double> z = ConvertComplex<double>(zz);
+  if (abs(z) == 0.) return 0;
   auto n = static_cast<double>(ni);
   auto one = std::complex<double> (1, 0);
-  return round((
+  auto lost_digits = round((
       abs(imag(z)) - log(2.) - n * ( real( log( z/n) + sqrt(one
       - pow2(z/n)) - log (one + sqrt (one
       - pow2(z/n)))
       )))/ log(10.));
+  return lost_digits;
 }
 
 template <typename FloatType>
@@ -77,6 +79,7 @@ int getNStar(int nmax, std::complex<FloatType> z, const int valid_digits) {
   if (nmax == 0) nmax = 1;
   int nstar = nmax;
   auto z_dp = ConvertComplex<double>(z);
+  if (std::abs(z_dp) == 0.) return  nstar;
   int forwardLoss = evalKapteynNumberOfLostSignificantDigits(nmax, z_dp);
   int increment = static_cast<int>(std::ceil(
       std::max(4* std::pow(std::abs( z_dp), 1/3.0), 5.0)

tests/CMakeLists.txt

@@ -6,6 +6,7 @@ include_directories(${GTEST_INCLUDE_DIRS})
 set(LIBS ${LIBS} ${GTEST_LIBRARIES})
 set(LIBS ${LIBS} pthread)
 
+add_compile_options(-D_GLIBCXX_DEBUG)
 # -- Output tests in directory
 
 add_executable("test_near_field"

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -84,7 +84,7 @@ TEST(an_test, mpmath_generated_input) {
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : an_test_30digits) {
     parse2_mpmath_data(min_abs_tol, data, x, m, n, an_mp, re_abs_tol, im_abs_tol);
-    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+    auto Nstop = nmie::LeRu_near_field_cutoff(m * x)+1;
 
     nmie::MultiLayerMie<nmie::FloatType> ml_mie;
     ml_mie.SetLayersSize({x});
@@ -111,7 +111,7 @@ TEST(bn_test, mpmath_generated_input) {
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : bn_test_30digits) {
     parse2_mpmath_data(min_abs_tol, data, x, m, n, bn_mp, re_abs_tol, im_abs_tol);
-    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+    auto Nstop = nmie::LeRu_near_field_cutoff(m * x)+1;
 
     nmie::MultiLayerMie<nmie::FloatType> ml_mie;
     ml_mie.SetLayersSize({x});
@@ -138,7 +138,7 @@ TEST(zeta_psizeta_test, mpmath_generated_input) {
   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 = nmie::LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_near_field_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+135), D3(Nstop+135),
         PsiZeta(Nstop+135), Psi(Nstop);
@@ -175,7 +175,7 @@ TEST(zeta_psizeta_test, mpmath_generated_input) {
 //  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 = nmie::LeRu_cutoff(z)+10000;
+//    auto Nstop = nmie::LeRu_near_field_cutoff(z)+10000;
 //    if (n > Nstop) continue;
 //    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop),
 //        PsiZeta(Nstop), zeta(Nstop);
@@ -200,7 +200,7 @@ TEST(psizeta_test, mpmath_generated_input) {
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : psi_mul_zeta_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, PsiZeta_mp, re_abs_tol, im_abs_tol);
-    auto Nstop = nmie::LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_near_field_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
     nmie::evalDownwardD1(z, D1dr);
@@ -225,7 +225,7 @@ TEST(psi_test, mpmath_generated_input) {
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : psi_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, Psi_mp, re_abs_tol, im_abs_tol);
-    auto Nstop = nmie::LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_near_field_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+35), Psi(Nstop);
     nmie::evalDownwardD1(z, D1dr);
@@ -247,7 +247,7 @@ TEST(D3test, mpmath_generated_input) {
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : D3_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, D3_mp, re_abs_tol, im_abs_tol);
-    auto Nstop = nmie::LeRu_cutoff(z)+35;
+    auto Nstop = nmie::LeRu_near_field_cutoff(z)+35;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
     nmie::evalDownwardD1(z, D1dr);
     nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
@@ -270,7 +270,7 @@ TEST(D3test, mpmath_generated_input) {
   for (const auto &data : D1_test_30digits) {
     parse2_mpmath_data(min_abs_tol, data, x, m, n, D1_mp, re_abs_tol, im_abs_tol);
     z = m*x;
-    auto Nstop = nmie::LeRu_cutoff(z)+1;
+    auto Nstop = nmie::LeRu_near_field_cutoff(z)+1;
 //    auto Nstop = n;
     std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop), r;
     int valid_digits = 10;

tests/test_bulk_sphere.cc

@@ -55,7 +55,7 @@ TEST(BulkSphere, HandlesInput) {
   for (const auto &data : parameters_and_results) {
     auto x = std::get<0>(data);
     auto m = std::get<1>(data);
-//    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+//    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
 
     nmie.SetLayersSize({x});
     nmie.SetLayersIndex({m});

tests/test_near_field.cc

@@ -18,7 +18,7 @@ TEST(RunFieldCalculationPolar, HandlesInput) {
   nmie.SetLayersSize({0.099});
   nmie.SetLayersIndex({ {0.75,0}});
   nmie.RunMieCalculation();
-//  nmie.RunFieldCalculationPolar(2);
+  nmie.RunFieldCalculationPolar(2);
 }
 //TEST(BulkSphere, HandlesInput) {
 //  nmie::MultiLayerMie<nmie::FloatType> nmie;
@@ -47,7 +47,7 @@ TEST(RunFieldCalculationPolar, HandlesInput) {
 //  for (const auto &data : parameters_and_results) {
 //    auto x = std::get<0>(data);
 //    auto m = std::get<1>(data);
-////    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+////    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
 //
 //    nmie.SetLayersSize({x});
 //    nmie.SetLayersIndex({m});