remove early cutoff for an_ and bn_

guiapp/src/components/PlotNearField.vue

@@ -17,7 +17,7 @@ import {
 } from 'vue'
 // import { flexRowTitleStyle } from 'components/config'
 import { plotlyChart } from 'src/store/plot-runtime/state'
-// import { Data, DataTitle } from 'plotly.js-dist-min'
+import { PlotData, /*, DataTitle*/ } from 'plotly.js-dist-min'
 import {nearFieldPlane} from 'src/store/simulation-setup/state';
 
 
@@ -58,20 +58,12 @@ export default defineComponent({
     const nearFieldEH = computed( ()=>$store.state.plotRuntime.nearFieldEH)
     watch([nearFieldEk, nearFieldHk, nearFieldEH], ()=>{
       nearFieldPlot.data.length = 0
-      if (crossSection.value == nearFieldPlane.Ek) {
-        nearFieldPlot.data.push({ z: $store.state.plotRuntime.nearFieldEk,
-          type: 'heatmap',  colorscale: 'Jet', })
-      }
-      if (crossSection.value == nearFieldPlane.Hk) {
-        nearFieldPlot.data.push({ z: $store.state.plotRuntime.nearFieldHk,
-          type: 'heatmap',  colorscale: 'Jet', })
-      }
-      if (crossSection.value == nearFieldPlane.EH) {
-        nearFieldPlot.data.push({ z: $store.state.plotRuntime.nearFieldEH,
-          type: 'heatmap',  colorscale: 'Jet', })
-      }
-      // if (crossSection.value == nearFieldPlane.Hk) $store.commit('plotRuntime/setNearFieldHk')
-      // if (crossSection.value == nearFieldPlane.EH) $store.commit('plotRuntime/setNearFieldEH')
+      const heatMapSettings: Partial<PlotData> = {type: 'heatmap',  colorscale: 'Jet'}
+      let heatMapData: Partial<PlotData> = {}
+      if (crossSection.value == nearFieldPlane.Ek) heatMapData  = { z: nearFieldEk.value}
+      if (crossSection.value == nearFieldPlane.Hk) heatMapData  = { z: nearFieldHk.value}
+      if (crossSection.value == nearFieldPlane.EH) heatMapData  = { z: nearFieldEH.value}
+      nearFieldPlot.data.push({...heatMapData, ...heatMapSettings})
 
     })
     // const xaxisTitle = computed(()=>{

guiapp/src/components/RunSimulationNearField.vue

@@ -118,6 +118,7 @@ export default defineComponent({
       }, 100)
     }
 
+    runNearFieldSimulation()
     watch(isRunning, ()=>{
           console.log(isRunning.value)
     })

guiapp/src/store/simulation-setup/state.ts

@@ -71,7 +71,7 @@ function setupFactory(hostIndex = 1,
                       plotLabel = '',
 
                       nearFieldSetup = {
-                        atWL: 619,
+                        atWL: 619.3885178,
                         relativePlotSize: 2,
                         plotSideResolution: 64,
                         crossSection: nearFieldPlane.Ek,

src/nmie-basic.hpp

@@ -459,8 +459,8 @@ namespace nmie {
                                                    std::vector<std::vector<FloatType> > &Pi,
                                                    std::vector<std::vector<FloatType> > &Tau) {
     const unsigned int perimeter_points = Pi.size();
-    for (auto &val:Pi) val.resize(available_maximal_nmax_);
-    for (auto &val:Tau) val.resize(available_maximal_nmax_);
+    for (auto &val:Pi) val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
+    for (auto &val:Tau) val.resize(available_maximal_nmax_, static_cast<FloatType>(0.0));
     double delta_Theta = eval_delta<double>(perimeter_points, from_Theta, to_Theta);
     for (unsigned int i=0; i < perimeter_points; i++) {
       auto Theta = static_cast<FloatType>(from_Theta + i*delta_Theta);
@@ -614,13 +614,13 @@ namespace nmie {
     std::vector<std::vector<std::complex<FloatType> > > Q(L), Ha(L), Hb(L);
 
     for (int l = 0; l < L; l++) {
-      Q[l].resize(nmax_ + 1);
-      Ha[l].resize(nmax_);
-      Hb[l].resize(nmax_);
+      Q[l].resize(nmax_ + 1, static_cast<FloatType>(0.0));
+      Ha[l].resize(nmax_, static_cast<FloatType>(0.0));
+      Hb[l].resize(nmax_, static_cast<FloatType>(0.0));
     }
 
-    an_.resize(nmax_);
-    bn_.resize(nmax_);
+    an_.resize(nmax_, static_cast<FloatType>(0.0));
+    bn_.resize(nmax_, static_cast<FloatType>(0.0));
 
     std::vector<std::complex<FloatType> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
 
@@ -733,17 +733,21 @@ namespace nmie {
         bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
       }
       if (n == 0) {a0 = cabs(an_[0]); b0 = cabs(bn_[0]);}
-      if (n == nmax_ - 1 && nmax_preset_ <= 0) {
+      if (n == nmax_ - 1 && nmax_preset_ <= 0
+          && (cabs(an_[n]) / a0 > convergence_threshold_ &&
+              cabs(bn_[n]) / b0 > convergence_threshold_)) {
         std::cout << "Failed to converge in Mie series for nmax="<<nmax_ << std::endl;
         std::cout << "convergence threshold: "<< convergence_threshold_ << std::endl;
         std::cout << "Mie series a[nmax]/a[1]:" << cabs(an_[n]) / a0 << " and b[nmax]/b[1]:" << cabs(bn_[n]) / b0 << std::endl;
 
       }
-      if (cabs(an_[n]) / a0 < convergence_threshold_ &&
-          cabs(bn_[n]) / b0 < convergence_threshold_) {
-        if (nmax_preset_ <= 0) nmax_ = n;
-        break;
-      }
+
+//  // TODO seems to provide not enough terms for near-field calclulation.
+//      if (cabs(an_[n]) / a0 < convergence_threshold_ &&
+//          cabs(bn_[n]) / b0 < convergence_threshold_) {
+//        if (nmax_preset_ <= 0) nmax_ = n;
+//        break;
+//      }
 
       if (nmm::isnan(an_[n].real()) || nmm::isnan(an_[n].imag()) ||
           nmm::isnan(bn_[n].real()) || nmm::isnan(bn_[n].imag())
@@ -816,7 +820,7 @@ namespace nmie {
     S1_.swap(tmp1);
     S2_ = S1_;
     // Precalculate cos(theta) - gives about 5% speed up.
-    std::vector<FloatType> costheta(theta_.size(), 0.0);
+    std::vector<FloatType> costheta(theta_.size(), static_cast<FloatType>(0.0));
     for (unsigned int t = 0; t < theta_.size(); t++) {
       costheta[t] = nmm::cos(theta_[t]);
     }

src/nmie-nearfield.hpp

@@ -93,10 +93,10 @@ namespace nmie {
     cln_.resize(L + 1);
     dln_.resize(L + 1);
     for (int l = 0; l <= L; l++) {
-      aln_[l].resize(nmax_);
-      bln_[l].resize(nmax_);
-      cln_[l].resize(nmax_);
-      dln_[l].resize(nmax_);
+      aln_[l].resize(nmax_, static_cast<FloatType>(0.0));
+      bln_[l].resize(nmax_, static_cast<FloatType>(0.0));
+      cln_[l].resize(nmax_, static_cast<FloatType>(0.0));
+      dln_[l].resize(nmax_, static_cast<FloatType>(0.0));
     }
 
     // Yang, paragraph under eq. A3
@@ -354,17 +354,17 @@ namespace nmie {
     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);
+//      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);
+      E[1] +=  Ex*cos_t(Phi)*cos_t(Theta);
+      E[2] += -Ex*sin_t(Phi);
       const auto Hy = Ex;
-      H[0] += Hy*sin(Theta)*sin(Phi);
-      H[1] += Hy*cos(Theta)*sin(Phi);
-      H[2] += Hy*cos(Phi);
+      H[0] += Hy*sin_t(Theta)*sin_t(Phi);
+      H[1] += Hy*cos_t(Theta)*sin_t(Phi);
+      H[2] += Hy*cos_t(Phi);
     }
 
     if( (!isConvergedE[0] || !isConvergedE[1] ||!isConvergedE[2] ||
@@ -545,9 +545,9 @@ void MultiLayerMie<FloatType>::calcRadialOnlyDependantFunctions(const double fro
 
     // Skip if not enough terms in Mie series (i.e. required near field nmax > available terms )
     if (near_field_nmax > available_maximal_nmax_)  near_field_nmax = available_maximal_nmax_;
-    Psi[j].resize(near_field_nmax + 1); D1n[j].resize(near_field_nmax + 1);
-    Zeta[j].resize(near_field_nmax + 1); D3n[j].resize(near_field_nmax + 1);
-    PsiZeta[j].resize(near_field_nmax + 1);
+    Psi[j].resize(near_field_nmax + 1, static_cast<FloatType>(0.0)); D1n[j].resize(near_field_nmax + 1, static_cast<FloatType>(0.0));
+    Zeta[j].resize(near_field_nmax + 1, static_cast<FloatType>(0.0)); D3n[j].resize(near_field_nmax + 1, static_cast<FloatType>(0.0));
+    PsiZeta[j].resize(near_field_nmax + 1, static_cast<FloatType>(0.0));
     std::complex<FloatType> ml;
     GetIndexAtRadius(Rho, ml);
     auto z = Rho*ml;

src/special-functions-impl.hpp

@@ -329,7 +329,7 @@ template <typename FloatType>
 void evalDownwardD1 (const std::complex<FloatType> z,
                      std::vector<std::complex<FloatType> >& D1) {
   int nmax = D1.size() - 1;
-  int valid_digits = 10;
+  int valid_digits = 16;
 #ifdef MULTI_PRECISION
   valid_digits += MULTI_PRECISION;
 #endif

tests/test_near_field.cc

@@ -8,7 +8,7 @@ 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);
-  double r = 2*nmie.PI_*100/619;
+  nmie::FloatType r = 2*nmie.PI_*100/619;
 //  double r = 1500;
   nmie.SetLayersSize({r/2, r});
   nmie.SetLayersIndex({ {4.0,0}, {4.0,0}});