bulk sphere tests e,g,h,i,k - passes

src/nmie-impl.hpp

@@ -275,6 +275,10 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
 
+  int LeRu_cutoff(std::complex<FloatType> z) {
+    auto x = std::abs(static_cast<std::complex<double>>(z));
+    return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
+  }
 
   // ********************************************************************** //
   // Calculate calcNstop - equation (17)                                    //
@@ -756,14 +760,17 @@ 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--) {
+//    for (int n = nmax_ - 2; n >= 0; n--) {
+    FloatType x2 = pow2(x.back());
+      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());
+        Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real())*2.0/x2;
         // 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());
+            + bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag())*2.0/x2;
+//        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())
             + ((n1 + n1 + 1.0) / (n1 * (n1 + 1.0))) * (an_[n] * std::conj(bn_[n])).real());
@@ -802,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/mpmath_input_arguments.py

@@ -10,8 +10,8 @@ complex_arguments = [
     # [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'],
-    # [0.055, [1.5, 1],    0.101491, 1.131687e-05, 'g'],
+    # [100,   [1.33,1e-5], 2.101321, 2.096594, 'e'],
+    [0.055, [1.5, 1],    0.101491, 1.131687e-05, 'g'],
     # [0.056, [1.5, 1],   0.1033467, 1.216311e-05, 'h'],
     # [100,   [1.5, 1],    2.097502, 1.283697, 'i'],
     # [1,     [10,  10],   2.532993, 2.049405, 'k'],

tests/mpmath_riccati_bessel.py

@@ -57,17 +57,23 @@ def D3(n,z):
 # bulk sphere
 # Ovidio
 def an(n, x, m):
-    print(f'D1 = {D1(n, m*x)}\n\
-    psi_n = {psi(n,x)}\n\
-    psi_nm1 = {psi(n-1,x)}\n\
-    ksi_n = {ksi(n,x)}\n\
-    ksi_nm1 = {ksi(n-1,x)}\n\
-    ')
+    # print(f'D1 = {D1(n, m*x)}\n\
+    # psi_n = {psi(n,x)}\n\
+    # psi_nm1 = {psi(n-1,x)}\n\
+    # ksi_n = {ksi(n,x)}\n\
+    # ksi_nm1 = {ksi(n-1,x)}\n\
+    # ')
     return (
             ( ( D1(n, m*x)/m + n/x )*psi(n,x) - psi(n-1,x) ) /
             ( ( D1(n, m*x)/m + n/x )*ksi(n,x) - ksi(n-1,x) )
     )
 def bn(n, x, m):
+    # print(f'D1 = {D1(n, m*x)}\n\
+    # psi_n = {psi(n,x)}\n\
+    # psi_nm1 = {psi(n-1,x)}\n\
+    # ksi_n = {ksi(n,x)}\n\
+    # ksi_nm1 = {ksi(n-1,x)}\n\
+    # ')
     return (
             ( ( D1(n, m*x)*m + n/x ) * psi(n,x) - psi(n-1,x) ) /
             ( ( D1(n, m*x)*m + n/x ) * ksi(n,x) - ksi(n-1,x) )
@@ -103,7 +109,7 @@ def Qany(x, m, nmax, output_dps, Qdiff):
         diff = Qdiff(n,x,m)
         Qany += diff
         Qnext = mp.nstr(Qany, output_dps)
-        Qlist.append(diff)
+        Qlist.append(Qany)
         i += 1
         if Qprev !=Qnext: i = 0
         Qprev = Qnext

tests/mpmath_special_functions_test_generator.py

@@ -212,6 +212,7 @@ def main():
     # sf_evals.run_test(mrb.ksi, 'zeta', is_only_x=True)
 
     # sf_evals.run_test(mrb.an, 'an', is_xm=True)
+    # sf_evals.run_test(mrb.bn, 'bn', is_xm=True)
 
     # sf_evals.run_test(mrb.psi, 'psi')
     # sf_evals.run_test(mrb.psi_div_ksi, 'psi_div_ksi')
@@ -233,12 +234,13 @@ def main():
         nmax = int(x + 4.05*x**(1./3.) + 2)+2+28
         print(f"\n ===== test case: {test_case} =====", flush=True)
         print(f"x={x}, m={m}, N={nmax} \nQsca_ref = {Qsca_ref}    \tQext_ref = {Qext_ref}", flush=True)
-        # Qext_mp = mrb.Qext(x,m,nmax, output_dps)
-        # Qsca_mp = mrb.Qsca(x,m,nmax, output_dps)
-        # print(f"Qsca_mp  = {mp.nstr(Qsca_mp[-1],output_dps)}    \tQext_mp  = {mp.nstr(Qext_mp[-1],output_dps)}", flush=True)
-        # print(mp.nstr(Qsca_mp,output_dps))
-        # print(mp.nstr(Qext_mp,output_dps))
+        Qext_mp = mrb.Qext(x,m,nmax, output_dps)
+        Qsca_mp = mrb.Qsca(x,m,nmax, output_dps)
+        print(f"Qsca_mp  = {mp.nstr(Qsca_mp[-1],output_dps)}    \tQext_mp  = {mp.nstr(Qext_mp[-1],output_dps)}", flush=True)
+        print(mp.nstr(Qsca_mp,output_dps))
+        print(mp.nstr(Qext_mp,output_dps))
+
         # n=1
-        # print(f'n={n}, x={x}, m={m}\nan[{n}]={mp.nstr(mrb.an(n,x,m), output_dps)}')
+        # print(f'n={n}, x={x}, m={m}\nbn[{n}]={mp.nstr(mrb.bn(n,x,m), output_dps)}')
 
 main()

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -42,10 +42,6 @@ std::vector< std::complex<double>>
               {-0.56047,0.11206e+01 }});
 
 
-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);
-}
 
 void parse_mpmath_data(const double min_abs_tol, const std::tuple< std::complex<double>, int, std::complex<double>, double, double > data,
                        std::complex<double> &z, int &n, std::complex<double> &func_mp,
@@ -80,15 +76,15 @@ 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;
   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 = LeRu_cutoff(m*x)+1;
+    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
 
     nmie::MultiLayerMie<nmie::FloatType> ml_mie;
     ml_mie.SetLayersSize({x});
@@ -107,6 +103,33 @@ TEST(an_test, mpmath_generated_input) {
   }
 }
 
+//TEST(bn_test, DISABLED_mpmath_generated_input) {
+TEST(bn_test, mpmath_generated_input) {
+  double min_abs_tol = 3e-14, x;
+  std::complex<double> m, bn_mp;
+  int n;
+  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;
+
+    nmie::MultiLayerMie<nmie::FloatType> ml_mie;
+    ml_mie.SetLayersSize({x});
+    ml_mie.SetLayersIndex({m});
+    ml_mie.SetMaxTerms(Nstop);
+    ml_mie.calcScattCoeffs();
+//    auto an = ml_mie.GetAn();
+    auto bn = ml_mie.GetBn();
+
+    if (n > bn.size()) continue;
+    if (n == 0) continue;
+    EXPECT_NEAR(std::real(bn[n-1]), std::real(bn_mp), re_abs_tol)
+              << "Db at n=" << n << " Nstop="<< Nstop<<" m="<<m<<" x="<<x;
+    EXPECT_NEAR(std::imag(bn[n-1]), std::imag(bn_mp), im_abs_tol)
+              << "Db at n=" << n << " Nstop="<< Nstop<<" m="<<m<<" x="<<x;
+  }
+}
+
 TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
 //TEST(zeta_psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 2e-10;
@@ -115,7 +138,7 @@ TEST(zeta_psizeta_test, DISABLED_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 = LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+135), D3(Nstop+135),
         PsiZeta(Nstop+135), Psi(Nstop);
@@ -152,7 +175,7 @@ TEST(zeta_test, DISABLED_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 = LeRu_cutoff(z)+10000;
+    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);
@@ -177,7 +200,7 @@ TEST(psizeta_test, DISABLED_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 = LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
     nmie::evalDownwardD1(z, D1dr);
@@ -202,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 = LeRu_cutoff(z)+10000;
+    auto Nstop = nmie::LeRu_cutoff(z)+10000;
     if (n > Nstop) continue;
     std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+35), Psi(Nstop);
     nmie::evalDownwardD1(z, D1dr);
@@ -224,7 +247,7 @@ TEST(D3test, DISABLED_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 = LeRu_cutoff(z)+35;
+    auto Nstop = nmie::LeRu_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);
@@ -247,7 +270,7 @@ TEST(D3test, DISABLED_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 = LeRu_cutoff(z)+1;
+    auto Nstop = nmie::LeRu_cutoff(z)+1;
     std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop+135), r;
     int valid_digits = 14;
     int nstar = nmie::getNStar(Nstop, z, valid_digits);

tests/test_bulk_sphere.cc

@@ -39,24 +39,28 @@ TEST(BulkSphere, HandlesInput) {
       parameters_and_results
       {
           // x, {Re(m), Im(m)}, Qext, Qsca, test_name
-//          {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'},
+          {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.101491, 1.131687e-05, 'g'},
-//          {0.056, {1.5, 1},   0.1033467, 1.216311e-05, 'h'},
-//          {100,   {1.5, 1},    2.097502, 1.283697, 'i'},
+          {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'},
+          {1,     {10,  10},   2.532993, 2.049405, 'k'},
 //          {100,   {10,  10,},  2.071124, 1.836785, 'l'},
 //          {10000, {10,  10},   2.005914, 1.795393, 'm'},
       };
   for (const auto &data : parameters_and_results) {
-    nmie.SetLayersSize({std::get<0>(data)});
-    nmie.SetLayersIndex({std::get<1>(data)});
-//    nmie.SetMaxTerms(150);
+    auto x = std::get<0>(data);
+    auto m = std::get<1>(data);
+    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+
+    nmie.SetLayersSize({x});
+    nmie.SetLayersIndex({m});
+    nmie.SetMaxTerms(Nstop);
     nmie.RunMieCalculation();
     double Qext = static_cast<double>(nmie.GetQext());
     double Qsca = static_cast<double>(nmie.GetQsca());

tests/test_spec_functions_data.hpp

@@ -1696,6 +1696,50 @@ an_test_30digits
 {100.0,{1.33,0.00001},190,{9.68143517313324082875587812679e-74,-3.55321407196406464604668533115e-69},9.7e-104,3.6e-99},
 };
 
+// x, complex(m), n, complex(f(n,z)), abs_err_real, abs_err_imag
+std::vector< std::tuple< nmie::FloatType, std::complex<nmie::FloatType>, int, std::complex<nmie::FloatType>, nmie::FloatType, nmie::FloatType > >
+bn_test_30digits
+= {
+{100.0,{1.33,0.00001},0,{0.985225810581282386106515279646,-0.115942975052557136687276065271},9.9e-31,1.2e-31},
+{100.0,{1.33,0.00001},1,{0.984198205480549450275215509099,0.121293322291954009772148249386},9.8e-31,1.2e-31},
+{100.0,{1.33,0.00001},2,{0.984111091723216958565637724715,-0.120455115559497447344633607914},9.8e-31,1.2e-31},
+{100.0,{1.33,0.00001},3,{0.988394155353531283859970050705,0.103170460307943493953358340537},9.9e-31,1.0e-31},
+{100.0,{1.33,0.00001},4,{0.981667453534137346692579191848,-0.129758755863461947283890344056},9.8e-31,1.3e-31},
+{100.0,{1.33,0.00001},5,{0.994393001262226790446889305308,0.0690616463467823781057223170999},9.9e-31,6.9e-32},
+{100.0,{1.33,0.00001},6,{0.978437542112137054090555573963,-0.141031795768905906171277851593},9.8e-31,1.4e-31},
+{100.0,{1.33,0.00001},7,{0.998928171031509403495922773509,0.017009070995538727841352739296},1.0e-30,1.7e-32},
+{100.0,{1.33,0.00001},8,{0.975504852376236081318591439551,-0.150432865515568326835206331805},9.8e-31,1.5e-31},
+{100.0,{1.33,0.00001},9,{0.996282201848085575689618507733,-0.0542687022438915322567305590184},1.0e-30,5.4e-32},
+{100.0,{1.33,0.00001},10,{0.974073321620675057800759094699,-0.154728677960347441147765309351},9.7e-31,1.5e-31},
+{100.0,{1.33,0.00001},11,{0.978237944059192548594231106649,-0.143326255920670845943686249146},9.8e-31,1.4e-31},
+{100.0,{1.33,0.00001},13,{0.935844619697458777560258716366,-0.243490008173120709350112588222},9.4e-31,2.4e-31},
+{100.0,{1.33,0.00001},14,{0.972331325754838412398557609546,-0.160146936582394111373556585325},9.7e-31,1.6e-31},
+{100.0,{1.33,0.00001},16,{0.961819306142619170862478976406,-0.188684975302727942310403913307},9.6e-31,1.9e-31},
+{100.0,{1.33,0.00001},18,{0.92859518928460200345277780158,-0.255629764591847580719924173528},9.3e-31,2.6e-31},
+{100.0,{1.33,0.00001},20,{0.849657444401319126749899621004,-0.356253019889860105519580422041},8.5e-31,3.6e-31},
+{100.0,{1.33,0.00001},22,{0.704236208513023092255978958968,-0.455570124927034131210396434143},7.0e-31,4.6e-31},
+{100.0,{1.33,0.00001},24,{0.501324251763404467667356357969,-0.499260921808715126231821237988},5.0e-31,5.0e-31},
+{100.0,{1.33,0.00001},27,{0.496604762530847106960828438922,-0.499038117595638483014877769869},5.0e-31,5.0e-31},
+{100.0,{1.33,0.00001},30,{0.12304253816148267349042194581,-0.326510066238273359692707119974},1.2e-31,3.3e-31},
+{100.0,{1.33,0.00001},34,{0.0497891171514503456820849654784,-0.215446225721925212102947512086},5.0e-32,2.2e-31},
+{100.0,{1.33,0.00001},38,{0.158370875419852237723715026027,0.364074618963252545967571248969},1.6e-31,3.6e-31},
+{100.0,{1.33,0.00001},42,{0.43712691436743175683325707469,0.494705346485833262243003552539},4.4e-31,4.9e-31},
+{100.0,{1.33,0.00001},47,{0.910854143460976674102618669457,0.282722158969578104293101897967},9.1e-31,2.8e-31},
+{100.0,{1.33,0.00001},52,{0.913786530245113986015944002211,-0.278857601617065918450227517047},9.1e-31,2.8e-31},
+{100.0,{1.33,0.00001},58,{0.0245164944062270740942982543122,-0.150320991662129562190444364757},2.5e-32,1.5e-31},
+{100.0,{1.33,0.00001},65,{0.880752944304931919955487086855,0.323172473965529844997584426499},8.8e-31,3.2e-31},
+{100.0,{1.33,0.00001},72,{0.418802868613230158697717538116,-0.492683390350917725808536465523},4.2e-31,4.9e-31},
+{100.0,{1.33,0.00001},80,{0.929238307111716727713442580138,-0.254491612243552408606177194977},9.3e-31,2.5e-31},
+{100.0,{1.33,0.00001},89,{0.784728369270289588227980415367,0.410501204235871868964872680977},7.8e-31,4.1e-31},
+{100.0,{1.33,0.00001},99,{0.0614997790835917145927623914716,0.239252646978288399998097825249},6.1e-32,2.4e-31},
+{100.0,{1.33,0.00001},111,{0.0000009291505307798752064389668148,0.000598869872117898826784602679261},9.3e-37,6.0e-34},
+{100.0,{1.33,0.00001},123,{3.0443025822485605231328209693e-13,2.7899384455209304186471850341e-10},3.0e-43,2.8e-40},
+{100.0,{1.33,0.00001},137,{5.21299608123152498131692259421e-24,-6.84402573284702321939804021489e-20},5.2e-54,6.8e-50},
+{100.0,{1.33,0.00001},153,{2.83618184239157249527152425776e-37,-5.53756343809554199199860075633e-33},2.8e-67,5.5e-63},
+{100.0,{1.33,0.00001},170,{3.2482367443426195596935606806e-53,-7.31964262759169298268294529878e-49},3.2e-83,7.3e-79},
+{100.0,{1.33,0.00001},190,{3.48336528648469774012102686213e-74,-8.50376657441951054200136719691e-70},3.5e-104,8.5e-100},
+};
+
 // complex(z), n, complex(f(n,z)), abs_err_real, abs_err_imag
 std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >
 psi_mul_zeta_test_16digits