update convergence evaluation for near-field computations

src/nmie-nearfield.hpp

@@ -175,13 +175,13 @@ namespace nmie {
         print_count++;
         std::cout<< std::setprecision(print_precision)
                  << "Warning: Potentially unstable calculation of aln[0]["
-                 << n << "] = "<< aln_[0][n] <<std::endl;
+                 << n << "] = "<< aln_[0][n] << " which is expected to be exact zero!"<<std::endl;
       }
       if (cabs(bln_[0][n]) > 1e-10  && print_count < 2)  {
         print_count++;
         std::cout<< std::setprecision(print_precision)
                  << "Warning: Potentially unstable calculation of bln[0]["
-                 << n << "] = "<< bln_[0][n] <<std::endl;
+                 << n << "] = "<< bln_[0][n] << " which is expected to be exact zero!" <<std::endl;
       }
       aln_[0][n] = 0.0;
       bln_[0][n] = 0.0;
@@ -280,6 +280,8 @@ namespace nmie {
 
     isConvergedE = {false, false, false}, isConvergedH = {false, false, false};
 //    evalType E0 = 0, H0=0;
+    std::vector< std::complex<evalType> > Ediff_prev = {{0.,0.},{0.,0.},{0.,0.}},
+        Hdiff_prev = {{0.,0.},{0.,0.},{0.,0.}};
     for (unsigned int n = 0; n < nmax; n++) {
       int n1 = n + 1;
       auto rn = static_cast<evalType>(n1);
@@ -298,7 +300,7 @@ namespace nmie {
       auto cln = ConvertComplex<evalType>(cln_[l][n]);
       auto dln = ConvertComplex<evalType>(dln_[l][n]);
       for (int i = 0; i < 3; i++) {
-        if (isConvergedE[i] && isConvergedH[i]) continue;
+        if (isConvergedE[i] && isConvergedH[i]) continue; // TODO is it safe?
         Ediff = En*(      cln*M1o1n[i] - c_i*dln*N1e1n[i]
                          + c_i*aln*N3e1n[i] -     bln*M3o1n[i]);
         Hdiff = En*(     -dln*M1e1n[i] - c_i*cln*N1o1n[i]
@@ -311,16 +313,19 @@ namespace nmie {
                     << " (of total nmax = "<<nmax<<")!!!"<<std::endl;
           break;
         }
-        if (n!=0) {
-          if (cabs(Ediff) == 0) isConvergedE[i] = true;
-          if (cabs(Hdiff) == 0) isConvergedH[i] = true;
-          if (cabs(E[i]) != 0.)
-            if (cabs(Ediff)/cabs(E[i]) < nearfield_convergence_threshold_) isConvergedE[i] = true;
-          if (cabs(H[i]) != 0.)
-            if (cabs(Hdiff)/cabs(H[i]) < nearfield_convergence_threshold_) isConvergedH[i] = true;
+        if (n>0) {
+          if (
+              (cabs(Ediff_prev[i]) <= cabs(E[i]) * nearfield_convergence_threshold_)
+                  &&  (cabs(Ediff) <= cabs(E[i]) * nearfield_convergence_threshold_)
+              ) isConvergedE[i] = true;
+          if (
+              (cabs(Hdiff_prev[i]) <= cabs(H[i]) * nearfield_convergence_threshold_)
+                  &&  (cabs(Hdiff) <= cabs(H[i]) * nearfield_convergence_threshold_)
+              ) isConvergedH[i] = true;
         }
-        if (isConvergedE[i]) Ediff = c_zero;
-        if (isConvergedH[i]) Hdiff = c_zero;
+        Ediff_prev[i] = Ediff;
+        Hdiff_prev[i] = Hdiff;
+
         if ((!isConvergedH[i] || !isConvergedE[i]) && n==nmax-1 && GetFieldConvergence()) {
           std::cout<<"Econv:"<<cabs(Ediff)/cabs(E[i])<<" Hconv:"<<cabs(Hdiff)/cabs(H[i])<<std::endl;
 
@@ -331,6 +336,9 @@ namespace nmie {
           H[i] += Hdiff;
           continue;
         }
+        if (n == 0) {
+
+        }
         if (n1 == mode_n_) {
           if (mode_type_ == Modes::kElectric || mode_type_ == Modes::kAll) {
             E[i] += En*( -c_i*dln*N1e1n[i]
@@ -359,8 +367,6 @@ namespace nmie {
 
     // Add the incident field
     if(l==L) {
-//      using nmm::sin_t;
-//      using nmm::cos_t;
       const auto z = Rho*cos_t(Theta);
       const auto Ex = std::complex<evalType>(cos_t(z), sin_t(z));
       E[0] +=  Ex*cos_t(Phi)*sin_t(Theta);

tests/test_near_field.cc

@@ -8,19 +8,45 @@ TEST(RunFieldCalculationCartesian, HandlesInput) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
 //  EXPECT_THROW(nmie.RunFieldCalculationPolar(0), std::invalid_argument);
 //  EXPECT_THROW(nmie.RunFieldCalculationPolar(1,1,10,5), std::invalid_argument);
-  nmie::FloatType r = 2*nmie.PI_*100/619;
+  nmie::FloatType total_r = 2*nmie.PI_*1000/532;
 //  double r = 1500;
-  nmie.SetLayersSize({r/2, r});
-  nmie.SetLayersIndex({ {4.0,0}, {4.0,0}});
+  nmie.SetLayersSize({total_r/2, total_r});
+  nmie.SetLayersIndex({ {1.330,0}, {1.33,0}});
   nmie.RunMieCalculation();
-  int nmax = 21;
+  double relative_max_distance = 1e-10;
+//  nmie.SetModeNmaxAndType(3,-1);
+//  int nmax = 21;
+  nmie.RunFieldCalculationCartesian(2, 5, relative_max_distance, nmie::Planes::kEk,
+                                    1.0, 0, 0, false,3);
+  auto Eabs = nmie.GetFieldEabs();
+  auto E = nmie.GetFieldE();
+  std::cout<<std::endl;
+  {
+    // Eabs points are located near the sphere outer border
+    //
+    //    0   1   2   3   4
+    //    ----- border ----
+    //    5   6   7   8   9
+    // distance between points (0) and (4) is relative_max_distance*total_r, initial
+    // value used for the test was 1e-10*total_r, so we expect good linear dependence
+    // for points from 0 to 4 and 5 to 9. In the asserts we check, that the slope doesn't
+    // change too fast inside the curve. While writing this, the test was failing.
+    // The value of z-coordinates of 2 and 7 points = 0
+    using nmie::nmm::abs;
+    EXPECT_TRUE(
+        ( abs(Eabs[0] - Eabs[1]) + abs(Eabs[3] - Eabs[4]) ) >= abs(Eabs[1] - Eabs[2])
+        );
+    EXPECT_TRUE(
+        ( abs(Eabs[5] - Eabs[6]) + abs(Eabs[8] - Eabs[9]) ) >= abs(Eabs[6] - Eabs[7])
+    );
+  }
+
+
+//  nmie.RunFieldCalculationCartesian(2, 2, 2, nmie::Planes::kHk,
+//                                    0, 0, 0, true);
+//  nmie.RunFieldCalculationCartesian(2, 2, 2, nmie::Planes::kEH,
+//                                    0, 0, 0, true);
   // TODO add check of E and H symmetry for X and Y axis inversion
-  nmie.RunFieldCalculationCartesian(2, 2, 2, nmie::Planes::kEk,
-                                    0, 0, 0, true, nmax);
-  nmie.RunFieldCalculationCartesian(2, 2, 2, nmie::Planes::kHk,
-                                    0, 0, 0, true, nmax);
-  nmie.RunFieldCalculationCartesian(2, 2, 2, nmie::Planes::kEH,
-                                    0, 0, 0, true, nmax);
 
 //  EXPECT_EQ(1, nmie.GetMaxTerms());
 //  EXPECT_FALSE(nmie.GetFieldConvergence());