update

src/nmie-applied.hpp

@@ -65,7 +65,9 @@ namespace nmie {
       } while (false);
     }
     // Set parameters in applied units 
-    void SetWavelength(FloatType wavelength) {wavelength_ = wavelength;};
+    void SetWavelength(FloatType wavelength) {
+      this->MultiLayerMie<FloatType>::MarkUncalculated();
+      wavelength_ = wavelength;};
     // It is possible to set only a multilayer target to run calculaitons.
     // For many runs it can be convenient to separate target and coating layers.
     // Per layer
@@ -87,6 +89,8 @@ namespace nmie {
     void SetFieldPointsSP(const std::vector< std::vector<FloatType> >& coords_sp);
 
     // Set common parameters
+    void SetModeNmaxAndType(int mode_n, int mode_type) {
+        this->MultiLayerMie<FloatType>::SetModeNmaxAndType(mode_n, mode_type);};
     void SetAnglesForPattern(FloatType from_angle, FloatType to_angle, int samples);
     std::vector<FloatType> GetAngles();
     

src/nmie-impl.cc

@@ -53,6 +53,7 @@
 #include <algorithm>
 #include <cstdio>
 #include <cstdlib>
+#include <cmath>
 #include <iostream>
 #include <iomanip>
 #include <stdexcept>
@@ -768,10 +769,10 @@ namespace nmie {
 
     const std::vector<FloatType>& x = size_param_;
 
-    MarkUncalculated();
+    //MarkUncalculated();
 
     // Calculate scattering coefficients
-    calcScattCoeffs();
+    if (!isScaCoeffsCalc_) calcScattCoeffs();
 
     // Initialize the scattering parameters
     Qext_ = 0.0;
@@ -787,6 +788,11 @@ namespace nmie {
     std::vector<std::complex<FloatType> > tmp1(theta_.size(),std::complex<FloatType>(0.0, 0.0));
     S1_.swap(tmp1);
     S2_ = S1_;
+    // Precalculate cos(theta) - gives about 5% speed up.
+    std::vector<FloatType> costheta(theta_.size(), 0.0);
+    for (unsigned int t = 0; t < theta_.size(); t++) {
+      costheta[t] = nmm::cos(theta_[t]);
+    }
 
     std::vector<FloatType> Pi(nmax_), Tau(nmax_);
 
@@ -795,31 +801,50 @@ 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 i = nmax_ - 2; i >= 0; i--) {
-      const int n = i + 1;
-      if (mode_n_ != Modes::kAll && n != mode_n_) continue;
-      // Equation (27)
-      Qext_ += (n + n + 1.0)*(an_[i].real() + bn_[i].real());
-      // Equation (28)
-      Qsca_ += (n + n + 1.0)*(an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
-                            + bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
-      // Equation (29)
-      Qpr_ += ((n*(n + 2.0)/(n + 1.0))*((an_[i]*std::conj(an_[n]) + bn_[i]*std::conj(bn_[n])).real())
-               + ((n + n + 1.0)/(n*(n + 1.0)))*(an_[i]*std::conj(bn_[i])).real());
-      // Equation (33)
-      Qbktmp += (FloatType)(n + n + 1.0)*(1.0 - 2.0*(n % 2))*(an_[i]- bn_[i]);
-      // Calculate the scattering amplitudes (S1 and S2)    //
-      // Precalculate cos(theta) - gives about 5% speed up.
-      std::vector<FloatType> costheta(theta_.size(), 0.0);
-      for (unsigned int t = 0; t < theta_.size(); t++) {
-        costheta[t] = nmm::cos(theta_[t]);
+    for (int n = nmax_ - 2; n >= 0; n--) {
+      const int n1 = n + 1;
+      if (mode_n_ == Modes::kAll) {
+        // Equation (27)
+        Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real());
+        // 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());
+        // 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());
+        // Equation (33)
+        Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (an_[n] - bn_[n]);
+        // Calculate the scattering amplitudes (S1 and S2) Equations (25a) - (25b)
+        for (unsigned int t = 0; t < theta_.size(); t++) {
+          calcPiTau(costheta[t], Pi, Tau);
+          S1_[t] += calc_S1(n1, an_[n], bn_[n], Pi[n], Tau[n]);
+          S2_[t] += calc_S2(n1, an_[n], bn_[n], Pi[n], Tau[n]);
+        }
+        continue;
       }
-      // Equations (25a) - (25b)                            //
-      for (unsigned int t = 0; t < theta_.size(); t++) {
-        calcPiTau(costheta[t], Pi, Tau);
-
-        S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
-        S2_[t] += calc_S2(n, an_[i], bn_[i], Pi[i], Tau[i]);
+      if (n1 == mode_n_) {
+        if (mode_type_ == Modes::kElectric || mode_type_ == Modes::kAll) {
+          Qext_ += (n1 + n1 + 1.0) * (an_[n].real());
+          Qsca_ += (n1 + n1 + 1.0) * (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag());
+          Qpr_ += std::nan("");
+          Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (an_[n]);
+          for (unsigned int t = 0; t < theta_.size(); t++) {
+            calcPiTau(costheta[t], Pi, Tau);
+            S1_[t] += calc_S1(n1, an_[n], static_cast<std::complex<FloatType>>(0), Pi[n], Tau[n]);
+            S2_[t] += calc_S2(n1, an_[n], static_cast<std::complex<FloatType>>(0), Pi[n], Tau[n]);
+          }
+        }
+        if (mode_type_ == Modes::kMagnetic || mode_type_ == Modes::kAll) {
+          Qext_ += (n1 + n1 + 1.0) * (bn_[n].real());
+          Qsca_ += (n1 + n1 + 1.0) * (bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
+          Qpr_ += std::nan("");
+          Qbktmp += (FloatType) (n1 + n1 + 1.0) * (1.0 - 2.0 * (n1 % 2)) * (bn_[n]);
+          for (unsigned int t = 0; t < theta_.size(); t++) {
+            calcPiTau(costheta[t], Pi, Tau);
+            S1_[t] += calc_S1(n1, static_cast<std::complex<FloatType>>(0), bn_[n], Pi[n], Tau[n]);
+            S2_[t] += calc_S2(n1, static_cast<std::complex<FloatType>>(0), bn_[n], Pi[n], Tau[n]);
+          }
+        }
       }
     }
     FloatType x2 = pow2(x.back());

src/nmie-js-wrapper.cc

@@ -45,6 +45,7 @@ EMSCRIPTEN_BINDINGS (c) {
                 .constructor<>()
                 .function("SetWavelength", &nmie::MultiLayerMieApplied<double>::SetWavelength)
                 .function("AddTargetLayerReIm",&nmie::MultiLayerMieApplied<double>::AddTargetLayerReIm)
+                .function("SetModeNmaxAndType",&nmie::MultiLayerMieApplied<double>::SetModeNmaxAndType)
                 .function("ClearTarget",&nmie::MultiLayerMieApplied<double>::ClearTarget)
                 .function("RunMieCalculation",&nmie::MultiLayerMieApplied<double>::RunMieCalculation)
                 .function("GetQsca",&nmie::MultiLayerMieApplied<double>::GetQsca)

web/index.html

@@ -9,21 +9,25 @@
 </head>
 <body>
 WL from <input type="number" id="fromWL" step="any" value="400.0">
-to <input type="number" id="toWL" step="any" value="800">
-step <input type="number" id="stepWL" step="any" value="0.1"><br>
+to <input type="number" id="toWL" step="any" value="1000">
+step <input type="number" id="stepWL" step="any" value="5"><br>
 R: <input type="number" id="R" step="any" value="100">
 Re(n): <input type="number" id="reN" step="any" value="4">
-Im(n): <input type="number" id="imN" step="any" value="0"><br>
-
+Im(n): <input type="number" id="imN" step="any" value="0.01"><br>
 <button id="Evaluate" autofocus>Run Mie!</button> <span id="time"></span>
 <div id="PlotQsca" style="width:640px;height:580px;" ></div>
   <script src="/nmie.js"></script>
   <script>
     const range = (start, stop, step = 1) =>
-            Array(Math.ceil((stop - start) / step)).fill(start).map((x, y) => x + y * step)
+            Array(Math.ceil((stop - start) / step)).fill(start).map((x, y) => x + y * step);
+    if (!Array.prototype.last){
+      Array.prototype.last = function(){
+        return this[this.length - 1];
+      };
+    };
     Module.addOnPostRun(() => {
-      const nmie = new Module.nmie()
       Evaluate.onclick = () => {
+        const nmie = new Module.nmie()
         var t0 = performance.now();
         var fromWL = parseFloat(document.getElementById("fromWL").value);
         var toWL = parseFloat(document.getElementById("toWL").value);
@@ -35,13 +39,34 @@ Im(n): <input type="number" id="imN" step="any" value="0"><br>
         nmie.AddTargetLayerReIm(R, reN, imN)
         var WLs = range(fromWL,toWL,stepWL)
         var Qsca = [], Qabs = [];
+        var Qsca_n = [[],[]], Qabs_n = [[],[]];
+        var mode_n = range(1,4)
+        var mode_n1 = range(1,5)
+        var mode_types = range(0,2)
+        mode_types.forEach(function(mode_type) {
+          mode_n1.forEach(function (n) {
+            Qsca_n[mode_type].push([]);
+            Qabs_n[mode_type].push([]);
+          });
+        });
         WLs.forEach(function(WL) {
+          console.log(WL)
+          nmie.SetModeNmaxAndType(-1, -1);
           nmie.SetWavelength(WL);
           nmie.RunMieCalculation();
           Qsca.push(nmie.GetQsca());
           Qabs.push(nmie.GetQabs());
+          mode_n.forEach(function (n) {
+            mode_types.forEach(function(mode_type) {
+              console.log(mode_type,n-1,Qsca_n[mode_type,0,n-1])
+              nmie.SetModeNmaxAndType(n,mode_type);
+              nmie.RunMieCalculation();
+              // Qsca_n[mode_type,n-1].push(nmie.GetQsca());
+              // Qabs_n[mode_type,n-1].push(nmie.GetQabs());
+            });
+          });
         });
-        console.log('Qsca = ', Qsca);
+        console.log('Qsca_n = ', Qsca_n);
         var trace1 = {
           x: WLs,
           y: Qsca,