fixes and additional NaN checks

Makefile

@@ -20,7 +20,7 @@ all:
 source:
 	$(PYTHON) setup.py sdist $(COMPILE) --dist-dir=../
 
-ext: $(SRCDIR)/nmie.cc $(SRCDIR)/nmie-pybind11.cc $(SRCDIR)/pb11_wrapper.cc
+ext: $(SRCDIR)/nmie.cc $(SRCDIR)/nmie-pybind11.cc $(SRCDIR)/pb11_wrapper.cc $(CXX_NMIE_HEADERS)
 	$(PYTHON) setup.py build_ext --inplace
 
 install:

examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -36,7 +36,10 @@ from scattnlay import scattnlay
 import numpy as np
 from matplotlib import pyplot as plt
 
-npts = 151
+import inspect
+print("Using scattnlay from ", inspect.getfile(scattnlay))
+
+npts = 251
 factor = 3. # plot extent compared to sphere radius
 index_H2O = 1.33+0.j
 
@@ -44,18 +47,22 @@ WL = 0.532 #mkm
 total_r = 1 #mkm
 isMP = False
 # isMP = True
+
+# nmax = 230
+nmax = -1
+
 nm = 1.0 # pihost medium
-x = 2.0 * np.pi * np.array([total_r / 4.0 * 3.0, total_r], dtype=np.float64) / WL
-m = np.array((index_H2O, index_H2O), dtype=np.complex128) / nm
+x = 2.0 * np.pi * np.array([total_r], dtype=np.float64) / WL
+m = np.array((index_H2O), dtype=np.complex128) / nm
 
 print("x =", x)
 print("m =", m)
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
     np.array([x]), np.array([m]))
-print("Qsca = " + str(Qsca))
+print("Qsca = " + str(Qsca)+" terms = "+str(terms))
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
     np.array([x]), np.array([m]), mp=True)
-print("mp Qsca = " + str(Qsca))
+print("mp Qsca = " + str(Qsca)+" terms = "+str(terms))
 
 scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
 zero = np.zeros(npts*npts, dtype = np.float64)
@@ -68,20 +75,25 @@ coordY = zero
 terms, E, H = fieldnlay(
     np.array([x]), np.array([m]),
     coordX, coordY, coordZ,
-    mp=isMP
+    mp=isMP,
+    nmax=nmax
 )
 Ec = E[0, :, :]
 Er = np.absolute(Ec)
 Eabs2 = (Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
 Eabs_data = np.resize(Eabs2, (npts, npts))
 label = r'$|E|^2$'
-plt.imshow(Eabs_data,
-           cmap='jet',
+pos = plt.imshow(Eabs_data,
+           cmap='gnuplot',
+                 # cmap='jet',
            vmin=0., vmax=14
 
            )
-print(np.min(Eabs_data), np.max(Eabs_data))
+plt.colorbar(pos)
+print(np.min(Eabs_data), np.max(Eabs_data)," terms = "+str(terms))
 mp = ''
 if isMP: mp = '_mp'
-plt.savefig("R"+str(total_r)+"mkm"+mp+".jpg")
+plt.savefig("R"+str(total_r)+"mkm"+mp+".jpg",
+            # dpi=300
+            )
 # plt.show()

src/nmie-impl.hpp

@@ -286,6 +286,7 @@ namespace nmie {
   // ********************************************************************** //
   template <typename FloatType>
   void MultiLayerMie<FloatType>::calcNstop() {
+    //Wiscombe
     const FloatType& xL = size_param_.back();
     if (xL <= 8) {
       nmax_ = newround(xL + 4.0*pow(xL, 1.0/3.0) + 1);
@@ -294,6 +295,9 @@ namespace nmie {
     } else {
       nmax_ = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
     }
+    //Le Ru
+    auto Nstop = nmie::LeRu_cutoff(static_cast<double>(xL))+1;
+    if (Nstop > nmax_) nmax_ = Nstop;
   }
 
 
@@ -693,9 +697,10 @@ namespace nmie {
       if (nmm::isnan(an_[n].real()) || nmm::isnan(an_[n].imag()) ||
           nmm::isnan(bn_[n].real()) || nmm::isnan(bn_[n].imag())
           ) {
-        nmax_ = n;
         // TODO somehow notify Python users about it
-        std::cout << "nmax value was chaned due to unexpected error. New values is "<<nmax_<<std::endl;
+        std::cout << "nmax value was chaned due to unexpected error. New values is "<< n
+                  << " (was "<<nmax_<<")"<<std::endl;
+        nmax_ = n;
         break;
       }
 
@@ -1028,7 +1033,8 @@ namespace nmie {
     // Calculate angular functions Pi and Tau
     calcPiTau(nmm::cos(Theta), Pi, Tau);
 
-    for (int n = nmax_ - 2; n >= 0; n--) {
+//    for (int n = nmax_ - 2; n >= 0; n--) {
+    for (int n = 0; n < nmax_-1; n++) {
       int n1 = n + 1;
       auto rn = static_cast<FloatType>(n1);
 
@@ -1040,13 +1046,17 @@ namespace nmie {
       std::complex<FloatType> En = ipow[n1 % 4]
       *static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
       for (int i = 0; i < 3; i++) {
+        auto Ediff = En*(      cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
+                         + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
+        auto Hdiff = En*(     -dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
+                         + c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
+        if (nmm::isnan(Ediff.real()) || nmm::isnan(Ediff.imag()) ||
+            nmm::isnan(Hdiff.real()) || nmm::isnan(Hdiff.imag())
+            ) break;
         if (mode_n_ == Modes::kAll) {
           // electric field E [V m - 1] = EF*E0
-          E[i] += En*(      cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
-                      + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
-
-          H[i] += En*(     -dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
-                      + c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
+          E[i] += Ediff;
+          H[i] += Hdiff;
           continue;
         }
         if (n1 == mode_n_) {

src/nmie.cc

@@ -257,7 +257,7 @@ namespace nmie {
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
   int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
-    return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, -1, -1);
+    return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, -1, -1);
 
   }
   //**********************************************************************************//
@@ -413,7 +413,7 @@ namespace nmie {
       ml_mie.SetFieldCoords({ConvertVector<FloatType>(Xp_vec),
 	    ConvertVector<FloatType>(Yp_vec),
 	    ConvertVector<FloatType>(Zp_vec) });
-
+      if (nmax != -1) ml_mie.SetMaxTerms(nmax);
       ml_mie.SetModeNmaxAndType(mode_n, mode_type);
 
       ml_mie.RunFieldCalculation();

src/special-functions-impl.hpp

@@ -62,8 +62,9 @@ namespace nmie {
 // Note, that Kapteyn seems to be too optimistic (at least by 3 digits
 // in some cases) for forward recurrence, see D1test with WYang_data
 int evalKapteynNumberOfLostSignificantDigits(const int ni,
-                                             const std::complex<double> z) {
+                                             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);
   auto n = static_cast<double>(ni);
   auto one = std::complex<double> (1, 0);
   return round((
@@ -74,6 +75,7 @@ int evalKapteynNumberOfLostSignificantDigits(const int ni,
 }
 
 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);
   int forwardLoss = evalKapteynNumberOfLostSignificantDigits(nmax, z_dp);
@@ -320,17 +322,23 @@ void evalForwardD1 (const std::complex<FloatType> z,
 void evalDownwardD1 (const std::complex<FloatType> z,
                      std::vector<std::complex<FloatType> >& D1) {
   int nmax = D1.size() - 1;
+  int valid_digits = 10;
+#ifdef MULTI_PRECISION
+  valid_digits += MULTI_PRECISION;
+#endif
+  int nstar = nmie::getNStar(nmax, z, valid_digits);
+  D1.resize(nstar+1);
   // Downward recurrence for D1 - equations (16a) and (16b)
-  D1[nmax] = std::complex<FloatType>(0.0, 0.0);
+  D1[nstar] = std::complex<FloatType>(0.0, 0.0);
   std::complex<FloatType> c_one(1.0, 0.0);
   const std::complex<FloatType> zinv = std::complex<FloatType>(1.0, 0.0)/z;
-  for (unsigned int n = nmax; n > 0; n--) {
+  for (unsigned int n = nstar; n > 0; n--) {
     D1[n - 1] = static_cast<FloatType>(n)*zinv - c_one/
         (D1[n] + static_cast<FloatType>(n)*zinv);
   }
   // Use D1[0] from upward recurrence
   D1[0] = complex_cot(z);
-
+  D1.resize(nmax+1);
 //  printf("D1[0] = (%16.15g, %16.15g) z=(%16.15g,%16.15g)\n", D1[0].real(),D1[0].imag(),
 //         z.real(),z.imag());
 }

tests/mpmath_input_arguments.py

@@ -15,7 +15,8 @@ complex_arguments = [
     # [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'],
-    [100,   [10,  10,],  2.071124, 1.836785, 'l'],
+    # [100,   [10,  10,],  2.071124, 1.836785, 'l'],
+    [11.8104986977059896, [1.33,0], 0,0, 'nanojet 1mkm']
     # [1, [mp.pi,  1],   0, 0, 'pi'],
     # [1, [mp.pi,  -1],   0, 0, 'pi'],
     # [1, [mp.pi,  mp.pi],   0, 0, 'pi'],

tests/mpmath_special_functions_test_generator.py

@@ -211,8 +211,8 @@ def main():
     # # In literature Zeta or Ksi denote the Riccati-Bessel function of third kind.
     # 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.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')

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -271,10 +271,12 @@ TEST(D3test, mpmath_generated_input) {
     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;
-    std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop+135), r;
-    int valid_digits = 14;
-    int nstar = nmie::getNStar(Nstop, z, valid_digits);
+//    auto Nstop = n;
+    std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop), r;
+    int valid_digits = 10;
+    int nstar = nmie::getNStar(n, z, valid_digits);
     r.resize(nstar);
+    Db.resize(nstar);
     nmie::evalBackwardR(z,r);
     nmie::convertRtoD1(z, r, Db);
     if (n > Db.size()) continue;
@@ -295,7 +297,7 @@ TEST(D3test, mpmath_generated_input) {
 
 //TEST(D1test, DISABLED_WYang_data){
 TEST(D1test, WYang_data){
-  double abs_tol = 1e-9;
+  double abs_tol = 4e-10;
   int test_loss_digits = std::round(15 - std::log10(1/abs_tol));
   int Nstop = 131;
   std::vector<std::complex<nmie::FloatType>> Df(Nstop), Db(Nstop),Dold(Nstop), r;
@@ -308,12 +310,6 @@ TEST(D1test, WYang_data){
   r.resize(Nstop+1);
   nmie::evalForwardR(z, r);
   nmie::convertRtoD1(z, r, Df);
-// eval backward recurrence
-  int valid_digits = 6;
-  int nstar = nmie::getNStar(Nstop, z, valid_digits);
-  r.resize(nstar);
-  nmie::evalBackwardR(z,r);
-  nmie::convertRtoD1(z, r, Db);
 
 for (int i = 0; i < Dtest_n.size(); i++) {
   int n = Dtest_n[i];
@@ -325,6 +321,14 @@ for (int i = 0; i < Dtest_n.size(); i++) {
     EXPECT_NEAR(std::imag(Df[n]), std::imag(Dtest_D1[i]),
                 abs_tol) << "f at n=" << n << " lost digits = " << forward_loss_digits;
   }
+// eval backward recurrence
+  int valid_digits = 6;
+  int nstar = nmie::getNStar(n, z, valid_digits);
+  r.resize(nstar);
+  Db.resize(nstar);
+  nmie::evalBackwardR(z,r);
+  nmie::convertRtoD1(z, r, Db);
+
   EXPECT_NEAR(std::real(Db[n]), std::real(Dtest_D1[i]),
               abs_tol) << "b at n=" << n;
   EXPECT_NEAR(std::imag(Db[n]), std::imag(Dtest_D1[i]),