Compute incident field directly instead of summing spherical wave exp…

…ansion for it

src/nmie-nearfield.hpp

@@ -164,6 +164,13 @@ namespace nmie {
       }  // end PEC condition
     }  // end of all l
 
+    // Remove the incident field coefficients. Incident field
+    // will be added as exp(ikz) after summation of the series.
+    for (int n = 0; n < nmax_; n++) {
+      cln_[L][n] = c_zero;
+      dln_[L][n] = c_zero;
+    }
+
     int print_precision = 16;
 #ifdef MULTI_PRECISION
     print_precision = MULTI_PRECISION;
@@ -350,6 +357,23 @@ namespace nmie {
           nmm::isnan(Hdiff_ft.real()) || nmm::isnan(Hdiff_ft.imag())
           ) break;
     }  // end of for all n
+
+    const unsigned L = refractive_index_.size();
+    // Add the incident field
+    if(l==L) {
+        using nmm::sin;
+        using nmm::cos;
+        const auto z = Rho*cos(Theta);
+        const auto Ex = std::complex<FloatType>(cos(z), sin(z));
+        E[0] +=  Ex*cos(Phi)*sin(Theta);
+        E[1] +=  Ex*cos(Phi)*cos(Theta);
+        E[2] += -Ex*sin(Phi);
+        const auto Hy = Ex;
+        H[0] += Hy*sin(Theta)*sin(Phi);
+        H[1] += Hy*cos(Theta)*sin(Phi);
+        H[2] += Hy*cos(Phi);
+    }
+
     if( (!isConvergedE[0] || !isConvergedE[1] ||!isConvergedE[2] ||
         !isConvergedH[0] || !isConvergedH[1] ||!isConvergedH[2] ) && GetFieldConvergence()) {
       std::cout << "Field evaluation failed to converge an nmax = "<< nmax << std::endl;