increased nmax, enabled all tests (without large size params)

src/nmie-impl.hpp

@@ -275,9 +275,10 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
 
-  int LeRu_cutoff(std::complex<FloatType> z) {
-    auto x = std::abs(static_cast<std::complex<double>>(z));
+  int LeRu_cutoff(std::complex<double> z) {
+    auto x = std::abs(z);
     return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
+//    return 10000;
   }
 
   // ********************************************************************** //
@@ -318,7 +319,7 @@ namespace nmie {
         riM1 = 0;
       nmax_ = std::max(nmax_, riM1);
     }
-    nmax_ += 15;  // Final nmax_ value
+    nmax_ += 25;  // Final nmax_ value
     // nmax_ *= nmax_;
     // printf("using nmax %i\n", nmax_);
   }
@@ -760,16 +761,15 @@ namespace nmie {
     // By using downward recurrence we avoid loss of precision due to float rounding errors
     // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
     //      http://en.wikipedia.org/wiki/Loss_of_significance
-//    for (int n = nmax_ - 2; n >= 0; n--) {
-    FloatType x2 = pow2(x.back());
-      for (int n = 0; n < nmax_; n++) {
+    for (int n = nmax_ - 2; n >= 0; n--) {
+//      for (int n = 0; n < nmax_; n++) {
       const int n1 = n + 1;
       if (mode_n_ == Modes::kAll) {
         // Equation (27)
-        Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real())*2.0/x2;
+        Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real());
         // Equation (28)
         Qsca_ += (n1 + n1 + 1.0) * (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag()
-            + bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag())*2.0/x2;
+            + bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
 //        std::cout<<"n ="<< n1 << " ext:"<<Qext_ <<" sca:"<<Qsca_<<std::endl;
         // Equation (29)
         Qpr_ += ((n1 * (n1 + 2.0) / (n1 + 1.0)) * ((an_[n] * std::conj(an_[n1]) + bn_[n] * std::conj(bn_[n1])).real())
@@ -809,9 +809,9 @@ namespace nmie {
         }
       }
     }
-//    FloatType x2 = pow2(x.back());
-//    Qext_ = 2.0*(Qext_)/x2;                                 // Equation (27)
-//    Qsca_ = 2.0*(Qsca_)/x2;                                 // Equation (28)
+    FloatType x2 = pow2(x.back());
+    Qext_ = 2.0*(Qext_)/x2;                                 // Equation (27)
+    Qsca_ = 2.0*(Qsca_)/x2;                                 // Equation (28)
     Qpr_ = Qext_ - 4.0*(Qpr_)/x2;                           // Equation (29)
     Qabs_ = Qext_ - Qsca_;                                  // Equation (30)
     albedo_ = Qsca_/Qext_;                                  // Equation (31)

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -76,8 +76,8 @@ void parse2_mpmath_data(const nmie::FloatType min_abs_tol,
 
 template<class T> inline T pow2(const T value) {return value*value;}
 
-TEST(an_test, DISABLED_mpmath_generated_input) {
-//TEST(an_test, mpmath_generated_input) {
+//TEST(an_test, DISABLED_mpmath_generated_input) {
+TEST(an_test, mpmath_generated_input) {
   double min_abs_tol = 3e-14, x;
   std::complex<double> m, an_mp;
   int n;
@@ -130,8 +130,8 @@ TEST(bn_test, mpmath_generated_input) {
   }
 }
 
-TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
-//TEST(zeta_psizeta_test, mpmath_generated_input) {
+//TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
+TEST(zeta_psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 2e-10;
   std::complex<double> z, zeta_mp;
   int n;
@@ -166,34 +166,34 @@ TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
   }
 }
 
-// Old way to evaluate Zeta
-TEST(zeta_test, DISABLED_mpmath_generated_input) {
-//TEST(zeta_test, mpmath_generated_input) {
-  double min_abs_tol = 2e-5;
-  std::complex<double> 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 = nmie::LeRu_cutoff(z)+10000;
-    if (n > Nstop) continue;
-    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop),
-        PsiZeta(Nstop), zeta(Nstop);
-    nmie::evalDownwardD1(z, D1dr);
-    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
-    nmie::evalUpwardZeta(z, D3, zeta);
-    if (std::isnan(std::real(zeta[n])) || std::isnan(std::imag(zeta[n]))) continue;
-
-    EXPECT_NEAR(std::real(zeta[n]), std::real(zeta_mp), re_abs_tol)
-              << "zeta[n] at n=" << n << " Nstop="<< Nstop<<" z="<<z;
-    EXPECT_NEAR(std::imag(zeta[n]), std::imag(zeta_mp), im_abs_tol)
-              << "zeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
-  }
-}
-
-
-TEST(psizeta_test, DISABLED_mpmath_generated_input) {
-//TEST(psizeta_test, mpmath_generated_input) {
+// // Old way to evaluate Zeta
+// TEST(zeta_test, DISABLED_mpmath_generated_input) {
+// //TEST(zeta_test, mpmath_generated_input) {
+//  double min_abs_tol = 2e-5;
+//  std::complex<double> 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 = nmie::LeRu_cutoff(z)+10000;
+//    if (n > Nstop) continue;
+//    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop),
+//        PsiZeta(Nstop), zeta(Nstop);
+//    nmie::evalDownwardD1(z, D1dr);
+//    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
+//    nmie::evalUpwardZeta(z, D3, zeta);
+//    if (std::isnan(std::real(zeta[n])) || std::isnan(std::imag(zeta[n]))) continue;
+//
+//    EXPECT_NEAR(std::real(zeta[n]), std::real(zeta_mp), re_abs_tol)
+//              << "zeta[n] at n=" << n << " Nstop="<< Nstop<<" z="<<z;
+//    EXPECT_NEAR(std::imag(zeta[n]), std::imag(zeta_mp), im_abs_tol)
+//              << "zeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
+//  }
+//}
+
+
+//TEST(psizeta_test, DISABLED_mpmath_generated_input) {
+TEST(psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 9e-11;
   std::complex<double> z, PsiZeta_mp;
   int n;
@@ -239,8 +239,8 @@ TEST(psi_test, mpmath_generated_input) {
 }
 
 
-TEST(D3test, DISABLED_mpmath_generated_input) {
-//TEST(D3test, mpmath_generated_input) {
+//TEST(D3test, DISABLED_mpmath_generated_input) {
+TEST(D3test, mpmath_generated_input) {
   double min_abs_tol = 2e-11;
   std::complex<double> z, D3_mp;
   int n;

tests/test_bulk_sphere.cc

@@ -4,8 +4,8 @@
 
 // TODO fails for MP with 100 digits.
 #ifndef MULTI_PRECISION
-TEST(BulkSphere, DISABLED_ArgPi) {
-//TEST(BulkSphere, ArgPi) {
+//TEST(BulkSphere, DISABLED_ArgPi) {
+TEST(BulkSphere, ArgPi) {
   std::vector<double> WLs{50, 80, 100,200, 400}; //nm
   double host_index = 2.;
   double core_radius = 100.; //nm
@@ -42,25 +42,25 @@ TEST(BulkSphere, HandlesInput) {
           {0.099, {0.75,0}, 7.417859e-06, 7.417859e-06, 'a'},
           {0.101, {0.75,0}, 8.033538e-06, 8.033538e-06, 'b'},
           {10,    {0.75,0},     2.232265, 2.232265, 'c'},
-//          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
           {100,   {1.33,1e-5}, 2.101321, 2.096594, 'e'},
-//          {10000, {1.33,1e-5}, 2.004089, 1.723857, 'f'},
           {0.055, {1.5, 1},    0.10149104, 1.131687e-05, 'g'},
           {0.056, {1.5, 1},   0.1033467, 1.216311e-05, 'h'},
           {100,   {1.5, 1},    2.097502, 1.283697, 'i'},
-//          {10000, {1.5, 1},    2.004368, 1.236574, 'j'},
           {1,     {10,  10},   2.532993, 2.049405, 'k'},
+          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
 //          {100,   {10,  10,},  2.071124, 1.836785, 'l'},
+//          {10000, {1.33,1e-5}, 2.004089, 1.723857, 'f'},
+//          {10000, {1.5, 1},    2.004368, 1.236574, 'j'},
 //          {10000, {10,  10},   2.005914, 1.795393, 'm'},
       };
   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_cutoff(m*x)+1;
 
     nmie.SetLayersSize({x});
     nmie.SetLayersIndex({m});
-    nmie.SetMaxTerms(Nstop);
+//    nmie.SetMaxTerms(Nstop);
     nmie.RunMieCalculation();
     double Qext = static_cast<double>(nmie.GetQext());
     double Qsca = static_cast<double>(nmie.GetQsca());