split web interface from nmie-applied (#74)

* split web interface from nmie-applied

* add windows

* add windows to actions

* more build tests

* test asser

* fix assert

.github/workflows/cmake.yml

@@ -11,7 +11,7 @@ env:
   BUILD_TYPE: Release
 
 jobs:
-  macosPython:
+  macos_Python:
     runs-on: macOS-latest
     steps:
       - uses: actions/checkout@v3
@@ -23,6 +23,18 @@ jobs:
         working-directory: ${{github.workspace}}
         run: tox run
 
+  # windows_Python:
+  #   runs-on: windows-latest
+  #   steps:
+  #     - uses: actions/checkout@v3
+
+  #     - name: Install tox
+  #       run: pip3 install tox
+
+  #     - name: Python initial test
+  #       working-directory: ${{github.workspace}}
+  #       run: tox run
+
   ubuntu_Python_wo_Boost:
     runs-on: ubuntu-latest
     steps:
@@ -68,6 +80,10 @@ jobs:
         working-directory: ${{github.workspace}}/build
         run: ctest -C ${{env.BUILD_TYPE}} --output-on-failure
 
+      - name: isBuilding without cmake
+        working-directory: ${{github.workspace}}/examples
+        run: ./go-cc-examples.sh
+
   ctest_wo_Boost:
     runs-on: ubuntu-latest
     steps:

examples/example-get-Mie.cc

@@ -1,49 +1,48 @@
 //**********************************************************************************//
-//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com> // Copyright
+//    (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
 //                                                                                  //
-//    This file is part of scattnlay                                                //
+//    This file is part of scattnlay //
 //                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
+//    This program is free software: you can redistribute it and/or modify // it
+//    under the terms of the GNU General Public License as published by // the
+//    Free Software Foundation, either version 3 of the License, or // (at your
+//    option) any later version.                                           //
 //                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
+//    This program is distributed in the hope that it will be useful, // but
+//    WITHOUT ANY WARRANTY; without even the implied warranty of //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the // GNU
+//    General Public License for more details. //
 //                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite the following reference:                                       //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//    The only additional remark is that we expect that all publications //
+//    describing work using this software, or all commercial products // using
+//    it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
+//        a multilayered sphere," Computer Physics Communications, // vol. 180,
+//        Nov. 2009, pp. 2348-2354.                                       //
 //                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//    You should have received a copy of the GNU General Public License // along
+//    with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 //   This program returns expansion coefficents of Mie series
 #include <complex>
 #include <cstdio>
 #include <string>
-#include "../src/nmie-applied.hpp"
 #include "../src/nmie-applied-impl.hpp"
 #include "../src/nmie-precision.hpp"
 #include "./read-spectra.h"
 
 // template<class T> inline T pow2(const T value) {return value*value;}
-int main(int argc, char *argv[]) {
-  using namespace nmie ;
+int main(int argc, char* argv[]) {
+  using namespace nmie;
   try {
     read_spectra::ReadSpectra Si_index, Ag_index;
     read_spectra::ReadSpectra plot_core_index_, plot_TiN_;
     std::string core_filename("Si-int.txt");
-    //std::string core_filename("Ag.txt");
-    //std::string TiN_filename("TiN.txt");
+    // std::string core_filename("Ag.txt");
+    // std::string TiN_filename("TiN.txt");
     std::string TiN_filename("Ag-int.txt");
-    //std::string TiN_filename("Si.txt");
+    // std::string TiN_filename("Si.txt");
     std::string shell_filename(core_filename);
 
     nmie::MultiLayerMieApplied<nmie::FloatType> multi_layer_mie;
@@ -51,105 +50,113 @@ int main(int argc, char *argv[]) {
     const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
     const std::complex<double> index_Si = std::sqrt(epsilon_Si);
     const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
-    double WL=500; //nm
-    double core_width = 5.27; //nm Si
-    double inner_width = 8.22; //nm Ag
-    double outer_width = 67.91; //nm  Si
+    double WL = 500;             // nm
+    double core_width = 5.27;    // nm Si
+    double inner_width = 8.22;   // nm Ag
+    double outer_width = 67.91;  // nm  Si
     bool isSiAgSi = true;
     double delta_width = 25.0;
-    //bool isSiAgSi = false;
+    // bool isSiAgSi = false;
     if (isSiAgSi) {
-      core_width = 5.27; //nm Si
-      inner_width = 8.22; //nm Ag
-      outer_width = 67.91; //nm  Si
+      core_width = 5.27;    // nm Si
+      inner_width = 8.22;   // nm Ag
+      outer_width = 67.91;  // nm  Si
       multi_layer_mie.AddTargetLayer(core_width, index_Si);
       multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
-      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+      multi_layer_mie.AddTargetLayer(outer_width + delta_width, index_Si);
     } else {
-      inner_width = 31.93; //nm Ag
-      outer_width = 4.06; //nm  Si
+      inner_width = 31.93;  // nm Ag
+      outer_width = 4.06;   // nm  Si
       multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
-      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+      multi_layer_mie.AddTargetLayer(outer_width + delta_width, index_Si);
     }
 
-    for (int dd = 0; dd<50; ++dd) {
+    for (int dd = 0; dd < 50; ++dd) {
       delta_width = dd;
-    FILE *fp;
-    std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat";
-    fp = fopen(fname.c_str(), "w");
+      FILE* fp;
+      std::string fname =
+          "absorb-layered-spectra-d" + std::to_string(dd) + ".dat";
+      fp = fopen(fname.c_str(), "w");
 
-    multi_layer_mie.SetWavelength(WL);
-    multi_layer_mie.RunMieCalculation();
-
-    double Qabs = static_cast<double>(multi_layer_mie.GetQabs());
-    printf("Qabs = %g\n", Qabs);
-    std::vector< std::vector<std::complex<nmie::FloatType> > > aln, bln, cln, dln;
-    multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
-    std::vector< std::vector<std::complex<double> > > d_aln =
-      nmie::ConvertComplexVectorVector<double>(aln);
-    std::string str = std::string("#WL ");
-    for (int l = 0; l<d_aln.size(); ++l) {
-      for (int n = 0; n<3; ++n) {
-	str+="|a|^2+|d|^2_ln"+std::to_string(l)+std::to_string(n)+" "
-	  + "|b|^2+|c|^2_ln"+std::to_string(l)+std::to_string(n)+" ";
-      }
-    }
-    str+="\n";
-    fprintf(fp, "%s", str.c_str());
-    fprintf(fp, "# |a|+|d|");
-    str.clear();
-
-    double from_WL = 400;
-    double to_WL = 600;
-    int total_points = 401;
-    double delta_WL = std::abs(to_WL - from_WL)/(total_points-1.0);
-    Si_index.ReadFromFile(core_filename).ResizeToComplex(from_WL, to_WL, total_points)
-      .ToIndex();
-    Ag_index.ReadFromFile(TiN_filename).ResizeToComplex(from_WL, to_WL, total_points)
-      .ToIndex();
-    auto Si_data = Si_index.GetIndex();
-    auto Ag_data = Ag_index.GetIndex();
-    for (int i=0; i < Si_data.size(); ++i) {
-      const double &WL = Si_data[i].first;
-      const std::complex<double> &Si = Si_data[i].second;
-      const std::complex<double> &Ag = Ag_data[i].second;
-      str+=std::to_string(WL);
-      multi_layer_mie.ClearTarget();
-      if (isSiAgSi) {
-	multi_layer_mie.AddTargetLayer(core_width, Si);
-	multi_layer_mie.AddTargetLayer(inner_width, Ag);
-	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
-      } else {
-	inner_width = 31.93; //nm Ag
-	outer_width = 4.06; //nm  Si
-	multi_layer_mie.AddTargetLayer(inner_width, Ag);
-	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
-      }
       multi_layer_mie.SetWavelength(WL);
       multi_layer_mie.RunMieCalculation();
-      multi_layer_mie.GetQabs();
-      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
-      for (int l = 0; l<aln.size(); ++l) {
-	for (int n = 0; n<3; ++n) {
-	  str+=" "+std::to_string(static_cast<double>(pow2(std::abs(aln[l][n]))+
-                                                      pow2(std::abs(dln[l][n]))))
-	    + " "
-	    + std::to_string(static_cast<double>(pow2(std::abs(bln[l][n]))
-                                                 + pow2(std::abs(cln[l][n])) ));
 
-	  // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n]))
-	  // 			  +dln[l][n].real() - pow2(std::abs(dln[l][n])))
-	  //   + " "
-	  //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
-	  // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n])) );
-	}
+      double Qabs = static_cast<double>(multi_layer_mie.GetQabs());
+      printf("Qabs = %g\n", Qabs);
+      std::vector<std::vector<std::complex<nmie::FloatType> > > aln, bln, cln,
+          dln;
+      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+      std::vector<std::vector<std::complex<double> > > d_aln =
+          nmie::ConvertComplexVectorVector<double>(aln);
+      std::string str = std::string("#WL ");
+      for (int l = 0; l < d_aln.size(); ++l) {
+        for (int n = 0; n < 3; ++n) {
+          str += "|a|^2+|d|^2_ln" + std::to_string(l) + std::to_string(n) +
+                 " " + "|b|^2+|c|^2_ln" + std::to_string(l) +
+                 std::to_string(n) + " ";
+        }
       }
-      str+="\n";
+      str += "\n";
       fprintf(fp, "%s", str.c_str());
+      fprintf(fp, "# |a|+|d|");
       str.clear();
-    }
 
-    fclose(fp);
+      double from_WL = 400;
+      double to_WL = 600;
+      int total_points = 401;
+      double delta_WL = std::abs(to_WL - from_WL) / (total_points - 1.0);
+      Si_index.ReadFromFile(core_filename)
+          .ResizeToComplex(from_WL, to_WL, total_points)
+          .ToIndex();
+      Ag_index.ReadFromFile(TiN_filename)
+          .ResizeToComplex(from_WL, to_WL, total_points)
+          .ToIndex();
+      auto Si_data = Si_index.GetIndex();
+      auto Ag_data = Ag_index.GetIndex();
+      for (int i = 0; i < Si_data.size(); ++i) {
+        const double& WL = Si_data[i].first;
+        const std::complex<double>& Si = Si_data[i].second;
+        const std::complex<double>& Ag = Ag_data[i].second;
+        str += std::to_string(WL);
+        multi_layer_mie.ClearTarget();
+        if (isSiAgSi) {
+          multi_layer_mie.AddTargetLayer(core_width, Si);
+          multi_layer_mie.AddTargetLayer(inner_width, Ag);
+          multi_layer_mie.AddTargetLayer(outer_width + delta_width, Si);
+        } else {
+          inner_width = 31.93;  // nm Ag
+          outer_width = 4.06;   // nm  Si
+          multi_layer_mie.AddTargetLayer(inner_width, Ag);
+          multi_layer_mie.AddTargetLayer(outer_width + delta_width, Si);
+        }
+        multi_layer_mie.SetWavelength(WL);
+        multi_layer_mie.RunMieCalculation();
+        multi_layer_mie.GetQabs();
+        multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+        for (int l = 0; l < aln.size(); ++l) {
+          for (int n = 0; n < 3; ++n) {
+            str += " " +
+                   std::to_string(static_cast<double>(
+                       pow2(std::abs(aln[l][n])) + pow2(std::abs(dln[l][n])))) +
+                   " " +
+                   std::to_string(static_cast<double>(
+                       pow2(std::abs(bln[l][n])) + pow2(std::abs(cln[l][n]))));
+
+            // str+=" "+std::to_string(aln[l][n].real() -
+            // pow2(std::abs(aln[l][n])) 			  +dln[l][n].real() -
+            // pow2(std::abs(dln[l][n])))
+            //   + " "
+            //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
+            // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n]))
+            // );
+          }
+        }
+        str += "\n";
+        fprintf(fp, "%s", str.c_str());
+        str.clear();
+      }
+
+      fclose(fp);
     }
 
     // WL = 500;
@@ -161,49 +168,57 @@ int main(int argc, char *argv[]) {
     // printf("\n Scattering");
     // for (int l = 0; l<aln.size(); ++l) {
     //   int n = 0;
-    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
-    //   n = 1;
-    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(),
+    //   aln[l][n].imag()); printf("bln[%i][%i] = %g, %gi\n", l, n+1,
+    //   bln[l][n].real(), bln[l][n].imag()); printf("cln[%i][%i] = %g, %gi\n",
+    //   l, n+1, cln[l][n].real(), cln[l][n].imag()); printf("dln[%i][%i] = %g,
+    //   %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag()); n = 1;
+    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(),
+    //   aln[l][n].imag()); printf("bln[%i][%i] = %g, %gi\n", l, n+1,
+    //   bln[l][n].real(), bln[l][n].imag()); printf("cln[%i][%i] = %g, %gi\n",
+    //   l, n+1, cln[l][n].real(), cln[l][n].imag()); printf("dln[%i][%i] = %g,
+    //   %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
     //   // n = 2;
-    //   // printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   // printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   // printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   // printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    //   // printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(),
+    //   aln[l][n].imag());
+    //   // printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(),
+    //   bln[l][n].imag());
+    //   // printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(),
+    //   cln[l][n].imag());
+    //   // printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(),
+    //   dln[l][n].imag());
     // }
     // printf("\n Absorbtion\n");
     // for (int l = 0; l<aln.size(); ++l) {
     //   if (l == aln.size()-1) printf(" Total ");
     //   printf("===== l=%i   =====\n", l);
     //   int n = 0;
-    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
-    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
-    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
-    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
-    //   n = 1;
-    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
-    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
-    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
-    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
+    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() -
+    //   pow2(std::abs(aln[l][n]))); printf("bln[%i][%i] = %g\n", l, n+1,
+    //   bln[l][n].real() - pow2(std::abs(bln[l][n]))); printf("cln[%i][%i] =
+    //   %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
+    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() -
+    //   pow2(std::abs(dln[l][n]))); n = 1; printf("aln[%i][%i] = %g\n", l, n+1,
+    //   aln[l][n].real() - pow2(std::abs(aln[l][n]))); printf("bln[%i][%i] =
+    //   %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
+    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() -
+    //   pow2(std::abs(cln[l][n]))); printf("dln[%i][%i] = %g\n", l, n+1,
+    //   dln[l][n].real() - pow2(std::abs(dln[l][n])));
     //   // n = 2;
-    //   // printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
-    //   // printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
-    //   // printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
-    //   // printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
+    //   // printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() -
+    //   pow2(std::abs(aln[l][n])));
+    //   // printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() -
+    //   pow2(std::abs(bln[l][n])));
+    //   // printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() -
+    //   pow2(std::abs(cln[l][n])));
+    //   // printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() -
+    //   pow2(std::abs(dln[l][n])));
     // }
 
-
-  } catch( const std::invalid_argument& ia ) {
+  } catch (const std::invalid_argument& ia) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
   }
-    return 0;
+  return 0;
 }
-
-

examples/example-minimal.cc

@@ -1,63 +1,69 @@
 //**********************************************************************************//
-//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com> // Copyright
+//    (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
 //                                                                                  //
-//    This file is part of scattnlay                                                //
+//    This file is part of scattnlay //
 //                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
+//    This program is free software: you can redistribute it and/or modify // it
+//    under the terms of the GNU General Public License as published by // the
+//    Free Software Foundation, either version 3 of the License, or // (at your
+//    option) any later version.                                           //
 //                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
+//    This program is distributed in the hope that it will be useful, // but
+//    WITHOUT ANY WARRANTY; without even the implied warranty of //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the // GNU
+//    General Public License for more details. //
 //                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite the following reference:                                       //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//    The only additional remark is that we expect that all publications //
+//    describing work using this software, or all commercial products // using
+//    it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
+//        a multilayered sphere," Computer Physics Communications, // vol. 180,
+//        Nov. 2009, pp. 2348-2354.                                       //
 //                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//    You should have received a copy of the GNU General Public License // along
+//    with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 //   This program evaluates absorption of a triple layered nanoparticle
+#include <cassert>
 #include <complex>
 #include <cstdio>
+#include <iomanip>
+#include <sstream>
 #include <string>
-#include "../src/nmie-applied.hpp"
+
 #include "../src/nmie-applied-impl.hpp"
 
-int main(int , char **) {
+int main(int, char**) {
   try {
     nmie::MultiLayerMieApplied<double> multi_layer_mie;
     const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
     const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
     const std::complex<double> index_Si = std::sqrt(epsilon_Si);
     const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
-    double WL=500; //nm
-    double core_width = 5.27; //nm Si
-    double inner_width = 8.22; //nm Ag
-    double outer_width = 67.91; //nm  Si
-    core_width = 5.27; //nm Si
-    inner_width = 8.22; //nm Ag
-    outer_width = 67.91; //nm  Si
+    double WL = 500;             // nm
+    double core_width = 5.27;    // nm Si
+    double inner_width = 8.22;   // nm Ag
+    double outer_width = 67.91;  // nm  Si
+    core_width = 5.27;           // nm Si
+    inner_width = 8.22;          // nm Ag
+    outer_width = 67.91;         // nm  Si
     multi_layer_mie.AddTargetLayer(core_width, index_Si);
     multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
     multi_layer_mie.AddTargetLayer(outer_width, index_Si);
     multi_layer_mie.SetWavelength(WL);
     multi_layer_mie.RunMieCalculation();
     double Qabs = multi_layer_mie.GetQabs();
-    printf("Qabs = %g\n", Qabs);
-  } catch( const std::invalid_argument &ia ) {
+    std::stringstream stream;
+    stream << std::fixed << std::setprecision(10) << Qabs;
+    auto Qabs_str = stream.str();
+    printf("Qabs = %s\n", Qabs_str.c_str());
+    assert(Qabs_str == "3.1415556911");
+
+  } catch (const std::invalid_argument& ia) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
   }
-    return 0;
+  return 0;
 }
-
-

examples/go-cc-examples.sh

@@ -10,23 +10,23 @@ PROGRAM='scattnlay-example.bin'
 # echo Compilation done. Running...
 # time ./$PROGRAM
 
-file=test-surf-integral.cc
-echo Compile $file with gcc
-rm -f $PROGRAM
-g++ -Ofast -std=c++11 $file  ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-echo Compilation done. Running...
-./$PROGRAM
-
-
-# file=example-minimal.cc
+# file=test-surf-integral.cc
 # echo Compile $file with gcc
 # rm -f $PROGRAM
-# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# g++ -Ofast -std=c++11 $file  ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# # g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
 # echo Compilation done. Running...
 # ./$PROGRAM
 
 
+file=example-minimal.cc
+echo Compile $file with gcc
+rm -f $PROGRAM
+g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc -lm  -o $PROGRAM -march=native -mtune=native -msse4.2
+echo Compilation done. Running...
+./$PROGRAM
+
+
 # file=example-get-Mie.cc
 # echo Compile $file with gcc
 # rm -f $PROGRAM

src/nmie-applied-impl.hpp

@@ -1,369 +1,389 @@
 #ifndef SRC_NMIE_APPLIED_IMPL_HPP_
 #define SRC_NMIE_APPLIED_IMPL_HPP_
-//**********************************************************************************//
-//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>          //
-//                                                                                  //
-//    This file is part of scattnlay                                                //
-//                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
-//                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
-//                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite at least one of the following references:                      //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
-//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie              //
-//        calculation of electromagnetic near-field for a multilayered              //
-//        sphere," Computer Physics Communications, vol. 214, May 2017,             //
-//        pp. 225-230.                                                              //
-//                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
-//                                                                                  //
-//    @brief  Wrapper class around nMie function for ease of use                    //
-//                                                                                  //
-//**********************************************************************************//
-#include <array>
+//**********************************************************************************
+//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite at least one of the following references:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+//        calculation of electromagnetic near-field for a multilayered
+//        sphere," Computer Physics Communications, vol. 214, May 2017,
+//        pp. 225-230.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+//    @brief  Wrapper class around nMie function for ease of use
+//
+//**********************************************************************************
 #include <algorithm>
+#include <array>
 #include <cstdio>
 #include <cstdlib>
 #include <stdexcept>
 #include <vector>
 
-namespace nmie {
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::GetFailed() {
-    FloatType faild_x = 9.42477796076938;
-    //FloatType faild_x = 9.42477796076937;
-    std::complex<FloatType> z(faild_x, 0.0);
-    std::vector<int> nmax_local_array = {20, 100, 500, 2500};
-    for (auto nmax_local : nmax_local_array) {
-      std::vector<std::complex<FloatType> > D1_failed(nmax_local + 1);
-      // Downward recurrence for D1 - equations (16a) and (16b)
-      D1_failed[nmax_local] = std::complex<FloatType>(0.0, 0.0);
-      const std::complex<FloatType> zinv = std::complex<FloatType>(1.0, 0.0)/z;
-      for (int n = nmax_local; n > 0; n--) {
-        D1_failed[n - 1] = FloatType(n)*zinv - 1.0/(D1_failed[n] + FloatType(n)*zinv);
-      }
-      printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n",
-             faild_x, nmax_local, D1_failed[0].real());
-    }
-    printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x,
-           calcD1confra(0,z).real());
-    //D1[nmax_] = calcD1confra(nmax_, z);
-
+#include "nmie-applied.hpp"
 
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::AddTargetLayer(FloatType width, std::complex<FloatType> layer_index) {
-      this->MarkUncalculated();
-      if (width <= 0)
-          throw std::invalid_argument("Layer width should be positive!");
-      target_width_.push_back(width);
-      target_index_.push_back(layer_index);
-  }  // end of void  MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::AddTargetLayerReIm(FloatType width,
-                                                           FloatType re_layer_index, FloatType im_layer_index) {
-    this->MarkUncalculated();
-    if (width <= 0)
-      throw std::invalid_argument("Layer width should be positive!");
-    target_width_.push_back(width);
-    target_index_.push_back(std::complex<FloatType>(re_layer_index,im_layer_index));
-  }  // end of void  MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetTargetPEC(FloatType radius) {
-    this->MarkUncalculated();
-    if (target_width_.size() != 0 || target_index_.size() != 0)
-      throw std::invalid_argument("Error! Define PEC target radius before any other layers!");
-    // Add layer of any index...
-    AddTargetLayer(radius, std::complex<FloatType>(0.0, 0.0));
-    // ... and mark it as PEC
-    this->SetPECLayer(0);
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<std::complex<FloatType> > index) {
-    this->MarkUncalculated();
-    coating_index_.clear();
-    for (auto value : index) coating_index_.push_back(value);
-  }  // end of void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<complex> index);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetCoatingWidth(std::vector<FloatType> width) {
-    this->MarkUncalculated();
-    coating_width_.clear();
-    for (auto w : width)
-      if (w <= 0)
-        throw std::invalid_argument("Coating width should be positive!");
-      else coating_width_.push_back(w);
-  }
-  // end of void MultiLayerMieApplied<FloatType>::SetCoatingWidth(...);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetWidthSP(const std::vector<FloatType>& size_parameter) {
-    this->MarkUncalculated();
-    this->size_param_.clear();
-    FloatType prev_size_parameter = 0.0;
-    for (auto layer_size_parameter : size_parameter) {
-      if (layer_size_parameter <= 0.0)
-        throw std::invalid_argument("Size parameter should be positive!");
-      if (prev_size_parameter > layer_size_parameter)
-        throw std::invalid_argument
-          ("Size parameter for next layer should be larger than the previous one!");
-      prev_size_parameter = layer_size_parameter;
-      this->size_param_.push_back(layer_size_parameter);
-    }
-  }
-  // end of void MultiLayerMieApplied<FloatType>::SetWidthSP(...);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetIndexSP(const std::vector< std::complex<FloatType> >& index) {
-    this->MarkUncalculated();
-    //refractive_index_.clear();
-    this->refractive_index_ = index;
-    // for (auto value : index) refractive_index_.push_back(value);
-  }  // end of void MultiLayerMieApplied<FloatType>::SetIndexSP(...);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(const std::vector< std::vector<FloatType> >& coords_sp) {
-    if (coords_sp.size() != 3)
-      throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
-    if (coords_sp[0].size() != coords_sp[1].size() || coords_sp[0].size() != coords_sp[2].size())
-      throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
-    this->coords_ = coords_sp;
-    // for (int i = 0; i < coords_sp_[0].size(); ++i) {
-    //   printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i], coords_sp_[2][i]);
-    // }
-  }  // end of void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::GenerateSizeParameter() {
-    this->MarkUncalculated();
-    this->size_param_.clear();
-    FloatType radius = 0.0;
-    for (auto width : target_width_) {
-      radius += width;
-      this->size_param_.push_back(2*nmie::PI_*radius/wavelength_);
-    }
-    for (auto width : coating_width_) {
-      radius += width;
-      this->size_param_.push_back(2*nmie::PI_*radius/wavelength_);
+namespace nmie {
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::GetFailed() {
+  FloatType faild_x = 9.42477796076938;
+  // FloatType faild_x = 9.42477796076937;
+  std::complex<FloatType> z(faild_x, 0.0);
+  std::vector<int> nmax_local_array = {20, 100, 500, 2500};
+  for (auto nmax_local : nmax_local_array) {
+    std::vector<std::complex<FloatType> > D1_failed(nmax_local + 1);
+    // Downward recurrence for D1 - equations (16a) and (16b)
+    D1_failed[nmax_local] = std::complex<FloatType>(0.0, 0.0);
+    const std::complex<FloatType> zinv = std::complex<FloatType>(1.0, 0.0) / z;
+    for (int n = nmax_local; n > 0; n--) {
+      D1_failed[n - 1] =
+          FloatType(n) * zinv - 1.0 / (D1_failed[n] + FloatType(n) * zinv);
     }
-    this->total_radius_ = radius;
-  }  // end of void MultiLayerMieApplied<FloatType>::GenerateSizeParameter();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::GenerateIndex() {
-    this->MarkUncalculated();
-    this->refractive_index_.clear();
-    for (auto index : this->target_index_)
-      this->refractive_index_.push_back(index);
-    for (auto index : this->coating_index_)
-      this->refractive_index_.push_back(index);
-  }  // end of void MultiLayerMieApplied<FloatType>::GenerateIndex();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius() {
-    if (!this->isMieCalculated())  GenerateSizeParameter();
-    return this->total_radius_;
-  }  // end of FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  std::vector< std::vector<FloatType> >
-  MultiLayerMieApplied<FloatType>::GetSpectra(FloatType from_WL, FloatType to_WL, int samples) {
-    if (!this->isMieCalculated())
-      throw std::invalid_argument("You should run calculations before result request!");
-    std::vector< std::vector<FloatType> > spectra;
-    FloatType step_WL = (to_WL - from_WL)/static_cast<FloatType>(samples);
-    FloatType wavelength_backup = wavelength_;
-    long fails = 0;
-    for (FloatType WL = from_WL; WL < to_WL; WL += step_WL) {
-      wavelength_ = WL;
-      try {
-        RunMieCalculation();
-      } catch(const std::invalid_argument& ia) {
-        fails++;
-        continue;
-      }
-      //printf("%3.1f ",WL);
-      spectra.push_back(std::vector<FloatType>({wavelength_, this->GetQext(),
-	      this->GetQsca(), this->GetQabs(), this->GetQbk()}));
-    }  // end of for each WL in spectra
-    printf("Spectrum has %li fails\n",fails);
-    wavelength_ = wavelength_backup;
-    return spectra;
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::ClearTarget() {
-    this->MarkUncalculated();
-    this->target_width_.clear();
-    this->target_index_.clear();
+    printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n",
+           faild_x, nmax_local, D1_failed[0].real());
   }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::ClearCoating() {
-    this->MarkUncalculated();
-    this->coating_width_.clear();
-    this->coating_index_.clear();
+  printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x,
+         calcD1confra(0, z).real());
+  // D1[nmax_] = calcD1confra(nmax_, z);
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::AddTargetLayer(
+    FloatType width,
+    std::complex<FloatType> layer_index) {
+  this->MarkUncalculated();
+  if (width <= 0)
+    throw std::invalid_argument("Layer width should be positive!");
+  target_width_.push_back(width);
+  target_index_.push_back(layer_index);
+}  // end of void  MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::AddTargetLayerReIm(
+    FloatType width,
+    FloatType re_layer_index,
+    FloatType im_layer_index) {
+  this->MarkUncalculated();
+  if (width <= 0)
+    throw std::invalid_argument("Layer width should be positive!");
+  target_width_.push_back(width);
+  target_index_.push_back(
+      std::complex<FloatType>(re_layer_index, im_layer_index));
+}  // end of void  MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetTargetPEC(FloatType radius) {
+  this->MarkUncalculated();
+  if (target_width_.size() != 0 || target_index_.size() != 0)
+    throw std::invalid_argument(
+        "Error! Define PEC target radius before any other layers!");
+  // Add layer of any index...
+  AddTargetLayer(radius, std::complex<FloatType>(0.0, 0.0));
+  // ... and mark it as PEC
+  this->SetPECLayer(0);
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetCoatingIndex(
+    std::vector<std::complex<FloatType> > index) {
+  this->MarkUncalculated();
+  coating_index_.clear();
+  for (auto value : index)
+    coating_index_.push_back(value);
+}  // end of void
+   // MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<complex>
+   // index);
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetCoatingWidth(
+    std::vector<FloatType> width) {
+  this->MarkUncalculated();
+  coating_width_.clear();
+  for (auto w : width)
+    if (w <= 0)
+      throw std::invalid_argument("Coating width should be positive!");
+    else
+      coating_width_.push_back(w);
+}
+// end of void MultiLayerMieApplied<FloatType>::SetCoatingWidth(...);
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetWidthSP(
+    const std::vector<FloatType>& size_parameter) {
+  this->MarkUncalculated();
+  this->size_param_.clear();
+  FloatType prev_size_parameter = 0.0;
+  for (auto layer_size_parameter : size_parameter) {
+    if (layer_size_parameter <= 0.0)
+      throw std::invalid_argument("Size parameter should be positive!");
+    if (prev_size_parameter > layer_size_parameter)
+      throw std::invalid_argument(
+          "Size parameter for next layer should be larger than the previous "
+          "one!");
+    prev_size_parameter = layer_size_parameter;
+    this->size_param_.push_back(layer_size_parameter);
   }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::ClearLayers() {
-    this->MarkUncalculated();
-    this->ClearTarget();
-    this->ClearCoating();
+}
+// end of void MultiLayerMieApplied<FloatType>::SetWidthSP(...);
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetIndexSP(
+    const std::vector<std::complex<FloatType> >& index) {
+  this->MarkUncalculated();
+  // refractive_index_.clear();
+  this->refractive_index_ = index;
+  // for (auto value : index) refractive_index_.push_back(value);
+}  // end of void MultiLayerMieApplied<FloatType>::SetIndexSP(...);
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(
+    const std::vector<std::vector<FloatType> >& coords_sp) {
+  if (coords_sp.size() != 3)
+    throw std::invalid_argument(
+        "Error! Wrong dimension of field monitor points!");
+  if (coords_sp[0].size() != coords_sp[1].size() ||
+      coords_sp[0].size() != coords_sp[2].size())
+    throw std::invalid_argument(
+        "Error! Missing coordinates for field monitor points!");
+  this->coords_ = coords_sp;
+  // for (int i = 0; i < coords_sp_[0].size(); ++i) {
+  //   printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i],
+  //   coords_sp_[2][i]);
+  // }
+}  // end of void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(...)
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::GenerateSizeParameter() {
+  this->MarkUncalculated();
+  this->size_param_.clear();
+  FloatType radius = 0.0;
+  for (auto width : target_width_) {
+    radius += width;
+    this->size_param_.push_back(2 * nmie::PI_ * radius / wavelength_);
   }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::ClearAllDesign() {
-    this->MarkUncalculated();
-    this->ClearLayers();
-    this->size_param_.clear();
-    this->refractive_index_.clear();
+  for (auto width : coating_width_) {
+    radius += width;
+    this->size_param_.push_back(2 * nmie::PI_ * radius / wavelength_);
   }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  //                         Computational core
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::ConvertToSP() {
-    this->MarkUncalculated();
-    if (target_width_.size() + coating_width_.size() == 0)
-      return;  // Nothing to convert, we suppose that SP was set directly
+  this->total_radius_ = radius;
+}  // end of void MultiLayerMieApplied<FloatType>::GenerateSizeParameter();
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::GenerateIndex() {
+  this->MarkUncalculated();
+  this->refractive_index_.clear();
+  for (auto index : this->target_index_)
+    this->refractive_index_.push_back(index);
+  for (auto index : this->coating_index_)
+    this->refractive_index_.push_back(index);
+}  // end of void MultiLayerMieApplied<FloatType>::GenerateIndex();
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius() {
+  if (!this->isMieCalculated())
     GenerateSizeParameter();
-    GenerateIndex();
-    if (this->size_param_.size() != this->refractive_index_.size())
-      throw std::invalid_argument("Ivalid conversion of width to size parameter units!/n");
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::RunMieCalculation() {
-    ConvertToSP();
-    this->MultiLayerMie<FloatType>::RunMieCalculation();
-  }
-
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::RunFieldCalculationPolar(const int outer_arc_points,
-                                                                 const int radius_points,
-                                                                 const double from_Rho, const double to_Rho,
-                                                                 const double from_Theta, const double to_Theta,
-                                                                 const double from_Phi, const double to_Phi,
-                                                                 const int isIgnoreAvailableNmax) {
-    ConvertToSP(); // Converts to size parameter units only the particle design,
-    // so we need to convert input parameters too...
-    const FloatType a = 2*nmie::PI_/wavelength_;
-    this->MultiLayerMie<FloatType>::RunFieldCalculationPolar(outer_arc_points, radius_points, a*from_Rho, a*to_Rho,
-                                                               from_Theta, to_Theta, from_Phi, to_Phi,
-                                                               isIgnoreAvailableNmax == 0 ? false : true);
-  }
-
+  return this->total_radius_;
+}  // end of FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius();
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+std::vector<std::vector<FloatType> >
+MultiLayerMieApplied<FloatType>::GetSpectra(FloatType from_WL,
+                                            FloatType to_WL,
+                                            int samples) {
+  if (!this->isMieCalculated())
+    throw std::invalid_argument(
+        "You should run calculations before result request!");
+  std::vector<std::vector<FloatType> > spectra;
+  FloatType step_WL = (to_WL - from_WL) / static_cast<FloatType>(samples);
+  FloatType wavelength_backup = wavelength_;
+  long fails = 0;
+  for (FloatType WL = from_WL; WL < to_WL; WL += step_WL) {
+    wavelength_ = WL;
+    try {
+      RunMieCalculation();
+    } catch (const std::invalid_argument& ia) {
+      fails++;
+      continue;
+    }
+    // printf("%3.1f ",WL);
+    spectra.push_back(
+        std::vector<FloatType>({wavelength_, this->GetQext(), this->GetQsca(),
+                                this->GetQabs(), this->GetQbk()}));
+  }  // end of for each WL in spectra
+  printf("Spectrum has %li fails\n", fails);
+  wavelength_ = wavelength_backup;
+  return spectra;
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
 template <typename FloatType>
-void MultiLayerMieApplied<FloatType>::RunFieldCalculationCartesian(const int first_side_points,
-                                                                   const int second_side_points,
-                                                            const double relative_side_length,
-                                                            const int plane_selected,
-                                                            const double at_x, const double at_y,
-                                                            const double at_z,
-                                                            const int isIgnoreAvailableNmax) {
-//  std::cout<<'test'<<std::endl;
+void MultiLayerMieApplied<FloatType>::ClearTarget() {
+  this->MarkUncalculated();
+  this->target_width_.clear();
+  this->target_index_.clear();
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::ClearCoating() {
+  this->MarkUncalculated();
+  this->coating_width_.clear();
+  this->coating_index_.clear();
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::ClearLayers() {
+  this->MarkUncalculated();
+  this->ClearTarget();
+  this->ClearCoating();
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::ClearAllDesign() {
+  this->MarkUncalculated();
+  this->ClearLayers();
+  this->size_param_.clear();
+  this->refractive_index_.clear();
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+//                         Computational core
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::ConvertToSP() {
+  this->MarkUncalculated();
+  if (target_width_.size() + coating_width_.size() == 0)
+    return;  // Nothing to convert, we suppose that SP was set directly
+  GenerateSizeParameter();
+  GenerateIndex();
+  if (this->size_param_.size() != this->refractive_index_.size())
+    throw std::invalid_argument(
+        "Ivalid conversion of width to size parameter units!/n");
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::RunMieCalculation() {
   ConvertToSP();
-  this->MultiLayerMie<FloatType>::RunFieldCalculationCartesian( first_side_points, second_side_points,
-                                                                relative_side_length, plane_selected,
-                                                                at_x, at_y, at_z,
-                                                                isIgnoreAvailableNmax == 0 ? false : true);
-
+  this->MultiLayerMie<FloatType>::RunMieCalculation();
 }
 
-//from https://toughengineer.github.io/demo/dsp/fft-perf/
-template <typename FloatType=double>
-emscripten::val toJSFloat64Array(const std::vector<double> &v) {
-  emscripten::val view{ emscripten::typed_memory_view(v.size(), v.data()) };  // make a view of local object
-
-  auto result = emscripten::val::global("Float64Array").new_(v.size());  // make a JS typed array
-  result.call<void>("set", view);  // set typed array values "on the JS side" using the memory view
-
-  return result;
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::RunFieldCalculationPolar(
+    const int outer_arc_points,
+    const int radius_points,
+    const double from_Rho,
+    const double to_Rho,
+    const double from_Theta,
+    const double to_Theta,
+    const double from_Phi,
+    const double to_Phi,
+    const int isIgnoreAvailableNmax) {
+  ConvertToSP();  // Converts to size parameter units only the particle design,
+  // so we need to convert input parameters too...
+  const FloatType a = 2 * nmie::PI_ / wavelength_;
+  this->MultiLayerMie<FloatType>::RunFieldCalculationPolar(
+      outer_arc_points, radius_points, a * from_Rho, a * to_Rho, from_Theta,
+      to_Theta, from_Phi, to_Phi, isIgnoreAvailableNmax == 0 ? false : true);
 }
 
 template <typename FloatType>
-emscripten::val MultiLayerMieApplied<FloatType>::GetFieldEabs() {
-  auto Eabs = this->MultiLayerMie<FloatType>::GetFieldEabs();
-  return toJSFloat64Array(Eabs);
+void MultiLayerMieApplied<FloatType>::RunFieldCalculationCartesian(
+    const int first_side_points,
+    const int second_side_points,
+    const double relative_side_length,
+    const int plane_selected,
+    const double at_x,
+    const double at_y,
+    const double at_z,
+    const int isIgnoreAvailableNmax) {
+  //  std::cout<<'test'<<std::endl;
+  ConvertToSP();
+  this->MultiLayerMie<FloatType>::RunFieldCalculationCartesian(
+      first_side_points, second_side_points, relative_side_length,
+      plane_selected, at_x, at_y, at_z,
+      isIgnoreAvailableNmax == 0 ? false : true);
 }
+
+// ********************************************************************** //
 // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  template <typename FloatType>
-  void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( std::vector< std::vector<std::complex<FloatType> > >& aln, std::vector< std::vector<std::complex<FloatType> > >& bln, std::vector< std::vector<std::complex<FloatType> > >& cln, std::vector< std::vector<std::complex<FloatType> > >& dln) {
-    ConvertToSP();  // Case of call before running full Mie calculation.
-    // Calculate scattering coefficients an_ and bn_
-    this->calcScattCoeffs();
-    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
-    this->calcExpanCoeffs();
-    aln = this->aln_;
-    bln = this->bln_;
-    cln = this->cln_;
-    dln = this->dln_;
+// ********************************************************************** //
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::GetExpanCoeffs(
+    std::vector<std::vector<std::complex<FloatType> > >& aln,
+    std::vector<std::vector<std::complex<FloatType> > >& bln,
+    std::vector<std::vector<std::complex<FloatType> > >& cln,
+    std::vector<std::vector<std::complex<FloatType> > >& dln) {
+  ConvertToSP();  // Case of call before running full Mie calculation.
+  // Calculate scattering coefficients an_ and bn_
+  this->calcScattCoeffs();
+  // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
+  this->calcExpanCoeffs();
+  aln = this->aln_;
+  bln = this->bln_;
+  cln = this->cln_;
+  dln = this->dln_;
 
-  }  // end of void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( ...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
+}  // end of void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( ...)
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
 
 }  // end of namespace nmie
 #endif  // SRC_NMIE_APPLIED_IMPL_HPP_

src/nmie-applied.cc

@@ -34,7 +34,6 @@
 #include <stdexcept>
 #include <vector>
 
-#include "nmie-applied.hpp"
 #include "nmie-applied-impl.hpp"
 
 namespace nmie {

src/nmie-applied.hpp

@@ -1,35 +1,35 @@
 #ifndef SRC_NMIE_APPLIED_HPP_
 #define SRC_NMIE_APPLIED_HPP_
-//**********************************************************************************//
-//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>          //
-//                                                                                  //
-//    This file is part of scattnlay                                                //
-//                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
-//                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
-//                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite at least one of the following references:                      //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
-//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie              //
-//        calculation of electromagnetic near-field for a multilayered              //
-//        sphere," Computer Physics Communications, vol. 214, May 2017,             //
-//        pp. 225-230.                                                              //
-//                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
-//**********************************************************************************//
+//**********************************************************************************
+//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite at least one of the following references:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+//        calculation of electromagnetic near-field for a multilayered
+//        sphere," Computer Physics Communications, vol. 214, May 2017,
+//        pp. 225-230.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//**********************************************************************************
 
 #include <array>
 #include <complex>
@@ -37,161 +37,232 @@
 #include <iostream>
 #include <vector>
 
-#include "nmie.hpp"
 #include "nmie-basic.hpp"
 #include "nmie-nearfield.hpp"
+#include "nmie.hpp"
 
-#include <emscripten/bind.h>
-#include <emscripten/val.h>
+namespace nmie {
 
+int nMieApplied(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);
+int nMieApplied(const unsigned int L,
+                std::vector<double>& x,
+                std::vector<std::complex<double> >& m,
+                const unsigned int nTheta,
+                std::vector<double>& Theta,
+                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);
+int nMieApplied(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,
+                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);
+int nMieApplied(const unsigned int L,
+                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);
 
-namespace nmie {
+template <typename FloatType = double>
+class MultiLayerMieApplied : public MultiLayerMie<FloatType> {
+  // Will throw for any error!
+ public:
+  void RunMieCalculation();
+  void RunFieldCalculationPolar(const int outer_arc_points,
+                                const int radius_points,
+                                const double from_Rho,
+                                const double to_Rho,
+                                const double from_Theta,
+                                const double to_Theta,
+                                const double from_Phi,
+                                const double to_Phi,
+                                const int isIgnoreAvailableNmax);
+  void RunFieldCalculationCartesian(const int first_side_points,
+                                    const int second_side_points,
+                                    const double relative_side_length,
+                                    const int plane_selected,
+                                    const double at_x,
+                                    const double at_y,
+                                    const double at_z,
+                                    const int isIgnoreAvailableNmax);
+  // emscripten::val GetFieldEabs();
+
+  void GetFailed();
+  long iformat = 0;
+  bool output = true;
+  void prn(FloatType var) {
+    do {
+      if (!output)
+        break;
+      ++iformat;
+      printf("%23.13e", var);
+      if (iformat % 4 == 0)
+        printf("\n");
+    } while (false);
+  }
+  // Set parameters in applied units
+  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
+  void AddTargetLayer(FloatType layer_width,
+                      std::complex<FloatType> layer_index);
+  void AddTargetLayerReIm(FloatType layer_width,
+                          FloatType re_layer_index,
+                          FloatType im_layer_index);
+  void AddCoatingLayer(FloatType layer_width,
+                       std::complex<FloatType> layer_index);
+  // For all layers
+  void SetTargetWidth(std::vector<FloatType> width);
+  void SetTargetIndex(std::vector<std::complex<FloatType> > index);
+  void SetTargetPEC(FloatType radius);
+  void SetCoatingWidth(std::vector<FloatType> width);
+  void SetCoatingIndex(std::vector<std::complex<FloatType> > index);
+  void SetFieldPoints(std::vector<std::array<FloatType, 3> > coords);
+
+  // Set parameters in size parameter units
+  void SetWidthSP(const std::vector<FloatType>& width);
+  void SetIndexSP(const std::vector<std::complex<FloatType> >& index);
+  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();
+
+  void ClearTarget();
+  void ClearCoating();
+  void ClearLayers();
+  void ClearAllDesign();  // Layers + SP + index_
+
+  // Applied units requests
+  FloatType GetTotalRadius();
+  FloatType GetTargetRadius();
+  FloatType GetCoatingWidth();
+  std::vector<FloatType> GetTargetLayersWidth();
+  std::vector<std::complex<FloatType> > GetTargetLayersIndex();
+  std::vector<FloatType> GetCoatingLayersWidth();
+  std::vector<std::complex<FloatType> > GetCoatingLayersIndex();
+  std::vector<std::vector<FloatType> > GetFieldPoints();
+  // ext, sca, abs, bk
+  std::vector<std::vector<FloatType> > GetSpectra(FloatType from_WL,
+                                                  FloatType to_WL,
+                                                  int samples);
+  FloatType GetRCSext();
+  FloatType GetRCSsca();
+  FloatType GetRCSabs();
+  FloatType GetRCSbk();
+  std::vector<FloatType> GetPatternEk();
+  std::vector<FloatType> GetPatternHk();
+  std::vector<FloatType> GetPatternUnpolarized();
+
+  // Dimensionless
+  FloatType GetQsca() { return this->MultiLayerMie<FloatType>::GetQsca(); };
+  FloatType GetQabs() { return this->MultiLayerMie<FloatType>::GetQabs(); };
+  FloatType GetQext() { return this->MultiLayerMie<FloatType>::GetQext(); };
+  // Size parameter units
+  std::vector<FloatType> GetLayerWidthSP();
+  // Same as to get target and coating index
+  std::vector<std::complex<FloatType> > GetLayerIndex();
+  std::vector<std::array<FloatType, 3> > GetFieldPointsSP();
+  // Do we need normalize field to size parameter?
+  /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldESP(); */
+  /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldHSP(); */
+  std::vector<std::array<FloatType, 5> > GetSpectraSP(
+      FloatType from_SP,
+      FloatType to_SP,
+      int samples);  // WL,ext, sca, abs, bk
+
+  std::vector<FloatType> GetPatternEkSP();
+  std::vector<FloatType> GetPatternHkSP();
+  std::vector<FloatType> GetPatternUnpolarizedSP();
+
+  void GetExpanCoeffs(std::vector<std::vector<std::complex<FloatType> > >& aln,
+                      std::vector<std::vector<std::complex<FloatType> > >& bln,
+                      std::vector<std::vector<std::complex<FloatType> > >& cln,
+                      std::vector<std::vector<std::complex<FloatType> > >& dln);
+
+  // Output results (data file + python script to plot it with matplotlib)
+  void PlotSpectra();
+  void PlotSpectraSP();
+  void PlotField();
+  void PlotFieldSP();
+  void PlotPattern();
+  void PlotPatternSP();
+
+ private:
+  void ConvertToSP();
+  void GenerateSizeParameter();
+  void GenerateIndex();
+  void InitMieCalculations();
+
+  void sbesjh(std::complex<FloatType> z,
+              std::vector<std::complex<FloatType> >& jn,
+              std::vector<std::complex<FloatType> >& jnp,
+              std::vector<std::complex<FloatType> >& h1n,
+              std::vector<std::complex<FloatType> >& h1np);
+  void sphericalBessel(std::complex<FloatType> z,
+                       std::vector<std::complex<FloatType> >& bj,
+                       std::vector<std::complex<FloatType> >& by,
+                       std::vector<std::complex<FloatType> >& bd);
+
+  FloatType wavelength_ = 1.0;
+  FloatType total_radius_ = 0.0;
+  /// Width and index for each layer of the structure
+  std::vector<FloatType> target_width_, coating_width_;
+  std::vector<std::complex<FloatType> > target_index_, coating_index_;
+
+  std::vector<std::vector<FloatType> > coords_sp_;
 
-  int nMieApplied(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);
-  int nMieApplied(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, 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);
-  int nMieApplied(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, 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);
-  int nMieApplied(const unsigned int L, 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);
-
-
-  template <typename FloatType = double>
-  class MultiLayerMieApplied : public MultiLayerMie<FloatType> {
-    // Will throw for any error!
-   public:
-    void RunMieCalculation();
-    void RunFieldCalculationPolar(const int outer_arc_points,
-                                  const int radius_points,
-                                  const double from_Rho, const double to_Rho,
-                                  const double from_Theta, const double to_Theta,
-                                  const double from_Phi, const double to_Phi,
-                                  const int isIgnoreAvailableNmax);
-    void RunFieldCalculationCartesian(const int first_side_points,
-                                      const int second_side_points,
-                                      const double relative_side_length,
-                                      const int plane_selected,
-                                      const double at_x, const double at_y,
-                                      const double at_z,
-                                      const int isIgnoreAvailableNmax);
-    emscripten::val GetFieldEabs();
-
-    void GetFailed();
-    long iformat = 0;
-    bool output = true;
-    void prn(FloatType var) {
-      do {
-	if (!output) break;
-	++iformat;
-	printf("%23.13e",var);
-	if (iformat%4 == 0) printf("\n");
-      } while (false);
-    }
-    // Set parameters in applied units
-    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
-    void AddTargetLayer(FloatType layer_width, std::complex<FloatType> layer_index);
-    void AddTargetLayerReIm(FloatType layer_width,
-            FloatType re_layer_index, FloatType im_layer_index);
-    void AddCoatingLayer(FloatType layer_width, std::complex<FloatType> layer_index);
-    // For all layers
-    void SetTargetWidth(std::vector<FloatType> width);
-    void SetTargetIndex(std::vector< std::complex<FloatType> > index);
-    void SetTargetPEC(FloatType radius);
-    void SetCoatingWidth(std::vector<FloatType> width);
-    void SetCoatingIndex(std::vector< std::complex<FloatType> > index);
-    void SetFieldPoints(std::vector< std::array<FloatType,3> > coords);
-
-    //Set parameters in size parameter units
-    void SetWidthSP(const std::vector<FloatType>& width);
-    void SetIndexSP(const std::vector< std::complex<FloatType> >& index);
-    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();
-
-    void ClearTarget();
-    void ClearCoating();
-    void ClearLayers();
-    void ClearAllDesign(); //Layers + SP + index_
-
-    // Applied units requests
-    FloatType GetTotalRadius();
-    FloatType GetTargetRadius();
-    FloatType GetCoatingWidth();
-    std::vector<FloatType>                  GetTargetLayersWidth();
-    std::vector< std::complex<FloatType> >  GetTargetLayersIndex();
-    std::vector<FloatType>                  GetCoatingLayersWidth();
-    std::vector< std::complex<FloatType> >  GetCoatingLayersIndex();
-    std::vector< std::vector<FloatType> >   GetFieldPoints();
-    std::vector< std::vector<FloatType> > GetSpectra(FloatType from_WL, FloatType to_WL, int samples);  // ext, sca, abs, bk
-    FloatType GetRCSext();
-    FloatType GetRCSsca();
-    FloatType GetRCSabs();
-    FloatType GetRCSbk();
-    std::vector<FloatType> GetPatternEk();
-    std::vector<FloatType> GetPatternHk();
-    std::vector<FloatType> GetPatternUnpolarized();
-
-    // Dimensionless
-    FloatType GetQsca(){ return this->MultiLayerMie<FloatType>::GetQsca();};
-    FloatType GetQabs(){ return this->MultiLayerMie<FloatType>::GetQabs();};
-    FloatType GetQext(){ return this->MultiLayerMie<FloatType>::GetQext();};
-    // Size parameter units
-    std::vector<FloatType> GetLayerWidthSP();
-    // Same as to get target and coating index
-    std::vector< std::complex<FloatType> > GetLayerIndex();
-    std::vector< std::array<FloatType,3> > GetFieldPointsSP();
-    // Do we need normalize field to size parameter?
-    /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldESP(); */
-    /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldHSP(); */
-    std::vector< std::array<FloatType,5> > GetSpectraSP(FloatType from_SP, FloatType to_SP, int samples);  // WL,ext, sca, abs, bk
-
-
-    std::vector<FloatType> GetPatternEkSP();
-    std::vector<FloatType> GetPatternHkSP();
-    std::vector<FloatType> GetPatternUnpolarizedSP();
-
-    void GetExpanCoeffs
-      (std::vector< std::vector<std::complex<FloatType> > >& aln,
-       std::vector< std::vector<std::complex<FloatType> > >& bln,
-       std::vector< std::vector<std::complex<FloatType> > >& cln,
-       std::vector< std::vector<std::complex<FloatType> > >& dln);
-
-
-    // Output results (data file + python script to plot it with matplotlib)
-    void PlotSpectra();
-    void PlotSpectraSP();
-    void PlotField();
-    void PlotFieldSP();
-    void PlotPattern();
-    void PlotPatternSP();
-
-  private:
-    void ConvertToSP();
-    void GenerateSizeParameter();
-    void GenerateIndex();
-    void InitMieCalculations();
-
-    void sbesjh(std::complex<FloatType> z, std::vector<std::complex<FloatType> >& jn,
-	            std::vector<std::complex<FloatType> >& jnp, std::vector<std::complex<FloatType> >& h1n,
-	            std::vector<std::complex<FloatType> >& h1np);
-    void sphericalBessel(std::complex<FloatType> z, std::vector<std::complex<FloatType> >& bj,
-			             std::vector<std::complex<FloatType> >& by, std::vector<std::complex<FloatType> >& bd);
-
-    FloatType wavelength_ = 1.0;
-    FloatType total_radius_ = 0.0;
-    /// Width and index for each layer of the structure
-    std::vector<FloatType> target_width_, coating_width_;
-    std::vector< std::complex<FloatType> > target_index_, coating_index_;
-
-    std::vector< std::vector<FloatType> > coords_sp_;
-
-  };  // end of class MultiLayerMie
+};  // end of class MultiLayerMie
 
 }  // end of namespace nmie
 #endif  // SRC_NMIE_APPLIED_HPP

src/nmie-js-wrapper.cc

@@ -1,100 +1,111 @@
-//**********************************************************************************//
-//    Copyright (C) 2009-2019  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2019  Konstantin Ladutenko <kostyfisik@gmail.com>          //
-//                                                                                  //
-//    This file is part of scattnlay                                                //
-//                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
-//                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
-//                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite at least one of the following references:                      //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
-//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie              //
-//        calculation of electromagnetic near-field for a multilayered              //
-//        sphere," Computer Physics Communications, vol. 214, May 2017,             //
-//        pp. 225-230.                                                              //
-//                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
-//                                                                                  //
+//**********************************************************************************
+//    Copyright (C) 2009-2019  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2019  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite at least one of the following references:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+//        calculation of electromagnetic near-field for a multilayered
+//        sphere," Computer Physics Communications, vol. 214, May 2017,
+//        pp. 225-230.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
 //    @brief  Wrapper to JS
-//                                                                                  //
-//**********************************************************************************//
-#include "nmie-applied.hpp"
-#include "nmie-applied-impl.hpp"
+//
+//**********************************************************************************
 #include "nmie-precision.hpp"
+#include "nmie-web-impl.hpp"
 
 using namespace emscripten;
 
-nmie::MultiLayerMieApplied<double> ml_mie;
+nmie::MultiLayerMieWeb<double> ml_mie;
 
-EMSCRIPTEN_BINDINGS (c) {
-    class_<nmie::MultiLayerMieApplied<double>>("nmie")
-        .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("RunFieldCalculationPolar",&nmie::MultiLayerMieApplied<double>::RunFieldCalculationPolar)
-        .function("RunFieldCalculationCartesian",&nmie::MultiLayerMieApplied<double>::RunFieldCalculationCartesian)
-        .function("GetFieldEabs",&nmie::MultiLayerMieApplied<double>::GetFieldEabs)
-        .function("GetQsca",&nmie::MultiLayerMieApplied<double>::GetQsca)
-        .function("GetQext",&nmie::MultiLayerMieApplied<double>::GetQext)
-        .function("GetQabs",&nmie::MultiLayerMieApplied<double>::GetQabs)
-//              .function("bf",&nmie::MultiLayerMieApplied<double>::bf)
-    ;
+EMSCRIPTEN_BINDINGS(c) {
+  class_<nmie::MultiLayerMieWeb<double>>("nmie")
+      .constructor<>()
+      .function("SetWavelength", &nmie::MultiLayerMieWeb<double>::SetWavelength)
+      .function("AddTargetLayerReIm",
+                &nmie::MultiLayerMieWeb<double>::AddTargetLayerReIm)
+      .function("SetModeNmaxAndType",
+                &nmie::MultiLayerMieWeb<double>::SetModeNmaxAndType)
+      .function("ClearTarget", &nmie::MultiLayerMieWeb<double>::ClearTarget)
+      .function("RunMieCalculation",
+                &nmie::MultiLayerMieWeb<double>::RunMieCalculation)
+      .function("RunFieldCalculationPolar",
+                &nmie::MultiLayerMieWeb<double>::RunFieldCalculationPolar)
+      .function("RunFieldCalculationCartesian",
+                &nmie::MultiLayerMieWeb<double>::RunFieldCalculationCartesian)
+      .function("GetFieldEabs", &nmie::MultiLayerMieWeb<double>::GetFieldEabs)
+      .function("GetQsca", &nmie::MultiLayerMieWeb<double>::GetQsca)
+      .function("GetQext", &nmie::MultiLayerMieWeb<double>::GetQext)
+      .function("GetQabs", &nmie::MultiLayerMieWeb<double>::GetQabs)
+      //              .function("bf",&nmie::MultiLayerMieWeb<double>::bf)
+      ;
 }
 
-//namespace nmie {
-  //**********************************************************************************//
-  // This function emulates a C call to calculate the actual scattering parameters    //
-  // and amplitudes.                                                                  //
-  //                                                                                  //
-  // Input parameters:                                                                //
-  //   L: Number of layers                                                            //
-  //   pl: Index of PEC layer. If there is none just send -1                          //
-  //   x: Array containing the size parameters of the layers [0..L-1]                 //
-  //   m: Array containing the relative refractive indexes of the layers [0..L-1]     //
-  //   nTheta: Number of scattering angles                                            //
-  //   Theta: Array containing all the scattering angles where the scattering         //
-  //          amplitudes will be calculated                                           //
-  //   nmax: Maximum number of multipolar expansion terms to be used for the          //
-  //         calculations. Only use it if you know what you are doing, otherwise      //
-  //         set this parameter to -1 and the function will calculate it              //
-  //                                                                                  //
-  // Output parameters:                                                               //
-  //   Qext: Efficiency factor for extinction                                         //
-  //   Qsca: Efficiency factor for scattering                                         //
-  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
-  //   Qbk: Efficiency factor for backscattering                                      //
-  //   Qpr: Efficiency factor for the radiation pressure                              //
-  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
-  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
-  //   S1, S2: Complex scattering amplitudes                                          //
-  //                                                                                  //
-  // Return value:                                                                    //
-  //   Number of multipolar expansion terms used for the calculations                 //
-  //**********************************************************************************//
-//  int nMieApplied(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) {
+// namespace nmie {
+//**********************************************************************************
+//  This function emulates a C call to calculate the actual scattering
+//  parameters and amplitudes.
+//
+//  Input parameters:
+//    L: Number of layers
+//    pl: Index of PEC layer. If there is none just send -1
+//    x: Array containing the size parameters of the layers [0..L-1]
+//    m: Array containing the relative refractive indexes of the layers [0..L-1]
+//    nTheta: Number of scattering angles
+//    Theta: Array containing all the scattering angles where the scattering
+//           amplitudes will be calculated
+//    nmax: Maximum number of multipolar expansion terms to be used for the
+//          calculations. Only use it if you know what you are doing, otherwise
+//          set this parameter to -1 and the function will calculate it
+//
+//  Output parameters:
+//    Qext: Efficiency factor for extinction
+//    Qsca: Efficiency factor for scattering
+//    Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)
+//    Qbk: Efficiency factor for backscattering
+//    Qpr: Efficiency factor for the radiation pressure
+//    g: Asymmetry factor (g = (Qext-Qpr)/Qsca)
+//    Albedo: Single scattering albedo (Albedo = Qsca/Qext)
+//    S1, S2: Complex scattering amplitudes
+//
+//  Return value:
+//    Number of multipolar expansion terms used for the calculations
+//**********************************************************************************
+//  int nMieWeb(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) {
 //
 //    if (x.size() != L || m.size() != L)
-//        throw std::invalid_argument("Declared number of layers do not fit x and m!");
+//        throw std::invalid_argument("Declared number of layers do not fit x
+//        and m!");
 //    if (Theta.size() != nTheta)
-//        throw std::invalid_argument("Declared number of sample for Theta is not correct!");
+//        throw std::invalid_argument("Declared number of sample for Theta is
+//        not correct!");
 //    try {
-//      MultiLayerMieApplied<FloatType> ml_mie;
+//      MultiLayerMieWeb<FloatType> ml_mie;
 //      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
 //      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
 //      ml_mie.SetAngles(ConvertVector<FloatType>(Theta));
@@ -122,17 +133,35 @@ EMSCRIPTEN_BINDINGS (c) {
 //    }
 //    return 0;
 //  }
-//  int nMieApplied(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, 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::nMieApplied(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+//  int nMieWeb(const unsigned int L, std::vector<double>& x,
+//  std::vector<std::complex<double> >& m, const unsigned int nTheta,
+//  std::vector<double>& Theta, 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::nMieWeb(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs,
+//    Qbk, Qpr, g, Albedo, S1, S2);
 //  }
-//  int nMieApplied(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, 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::nMieApplied(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+//  int nMieWeb(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, 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::nMieWeb(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs,
+//    Qbk, Qpr, g, Albedo, S1, S2);
 //  }
-//  int nMieApplied(const unsigned int L, 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::nMieApplied(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+//  int nMieWeb(const unsigned int L, 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::nMieWeb(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs,
+//    Qbk, Qpr, g, Albedo, S1, S2);
 //  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
 
 //}  // end of namespace nmie

src/nmie-precision.hpp

@@ -1,117 +1,130 @@
 #ifndef SRC_NMIE_PRECISION_H_
 #define SRC_NMIE_PRECISION_H_
 //**********************************************************************************//
-//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com> // Copyright
+//    (C) 2013-2018  Konstantin Ladutenko <kostyfisik@gmail.com>          //
 //                                                                                  //
-//    This file is part of scattnlay                                                //
+//    This file is part of scattnlay //
 //                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
+//    This program is free software: you can redistribute it and/or modify // it
+//    under the terms of the GNU General Public License as published by // the
+//    Free Software Foundation, either version 3 of the License, or // (at your
+//    option) any later version.                                           //
 //                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
+//    This program is distributed in the hope that it will be useful, // but
+//    WITHOUT ANY WARRANTY; without even the implied warranty of //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the // GNU
+//    General Public License for more details. //
 //                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite at least one of the following references:                      //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
-//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie              //
-//        calculation of electromagnetic near-field for a multilayered              //
-//        sphere," Computer Physics Communications, vol. 214, May 2017,             //
-//        pp. 225-230.                                                              //
+//    The only additional remark is that we expect that all publications //
+//    describing work using this software, or all commercial products // using
+//    it, cite at least one of the following references:                      //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
+//        a multilayered sphere," Computer Physics Communications, // vol. 180,
+//        Nov. 2009, pp. 2348-2354.                                       //
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie //
+//        calculation of electromagnetic near-field for a multilayered //
+//        sphere," Computer Physics Communications, vol. 214, May 2017, // pp.
+//        225-230. //
 //                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//    You should have received a copy of the GNU General Public License // along
+//    with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
 //**********************************************************************************//
-// This class implements the algorithm for a multilayered sphere described by:      //
-//    [1] W. Yang, "Improved recursive algorithm for light scattering by a          //
-//        multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp. 1710-1720.  //
+// This class implements the algorithm for a multilayered sphere described by:
+// //
+//    [1] W. Yang, "Improved recursive algorithm for light scattering by a //
+//        multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp.
+//        1710-1720.  //
 //                                                                                  //
-// You can find the description of all the used equations in:                       //
-//    [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+// You can find the description of all the used equations in: //
+//    [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
+//        a multilayered sphere," Computer Physics Communications, // vol. 180,
+//        Nov. 2009, pp. 2348-2354.                                       //
 //                                                                                  //
-// Hereinafter all equations numbers refer to [2]                                   //
+// Hereinafter all equations numbers refer to [2] //
 //**********************************************************************************//
+#include <complex>
+#include <vector>
 #ifdef MULTI_PRECISION
-#include <boost/multiprecision/number.hpp>
 #include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/number.hpp>
 #endif  // MULTI_PRECISION
 
 namespace nmie {
-  #ifdef MULTI_PRECISION
-  namespace nmm = boost::multiprecision;
-  typedef nmm::number<nmm::cpp_bin_float<MULTI_PRECISION> > FloatType;
-  #else
-  namespace nmm = std;
-  typedef double FloatType;
-  //typedef float FloatType;
-  #endif  // MULTI_PRECISION
-
-  template<class T> T sin_t(T v) {
-    if (std::is_same<T, double>::value) return static_cast<T>(std::sin(static_cast<double>(v)));
-    return static_cast<T>(nmm::sin(static_cast<FloatType >(v)));
-  }
-  template<class T> T cos_t(T v) {
-    if (std::is_same<T, double>::value) return static_cast<T>(std::cos(static_cast<double>(v)));
-    return static_cast<T>(nmm::cos(static_cast<FloatType >(v)));
-  }
-  template<class T> T sqrt_t(T v) {
-    if (std::is_same<T, double>::value) return static_cast<T>(std::sqrt(static_cast<double>(v)));
-    return static_cast<T>(nmm::sqrt(static_cast<FloatType >(v)));
-  }
-
+#ifdef MULTI_PRECISION
+namespace nmm = boost::multiprecision;
+typedef nmm::number<nmm::cpp_bin_float<MULTI_PRECISION> > FloatType;
+#else
+namespace nmm = std;
+typedef double FloatType;
+// typedef float FloatType;
+#endif  // MULTI_PRECISION
 
+template <class T>
+T sin_t(T v) {
+  if (std::is_same<T, double>::value)
+    return static_cast<T>(std::sin(static_cast<double>(v)));
+  return static_cast<T>(nmm::sin(static_cast<FloatType>(v)));
+}
+template <class T>
+T cos_t(T v) {
+  if (std::is_same<T, double>::value)
+    return static_cast<T>(std::cos(static_cast<double>(v)));
+  return static_cast<T>(nmm::cos(static_cast<FloatType>(v)));
+}
+template <class T>
+T sqrt_t(T v) {
+  if (std::is_same<T, double>::value)
+    return static_cast<T>(std::sqrt(static_cast<double>(v)));
+  return static_cast<T>(nmm::sqrt(static_cast<FloatType>(v)));
+}
 
 template <typename ToFloatType, typename FromFloatType>
-  std::vector<ToFloatType> ConvertVector(const std::vector<FromFloatType> x) {
-    std::vector<ToFloatType> new_x;
-    for (auto element : x) {
-      new_x.push_back(static_cast<ToFloatType>(element));
-    }
-    return new_x;
+std::vector<ToFloatType> ConvertVector(const std::vector<FromFloatType> x) {
+  std::vector<ToFloatType> new_x;
+  for (auto element : x) {
+    new_x.push_back(static_cast<ToFloatType>(element));
   }
+  return new_x;
+}
 
-  template <typename ToFloatType, typename FromFloatType>
-  std::complex<ToFloatType> ConvertComplex(std::complex<FromFloatType> z) {
-    return std::complex<ToFloatType>(static_cast<ToFloatType>(z.real()),
-                                     static_cast<ToFloatType>(z.imag()));
-  }
+template <typename ToFloatType, typename FromFloatType>
+std::complex<ToFloatType> ConvertComplex(std::complex<FromFloatType> z) {
+  return std::complex<ToFloatType>(static_cast<ToFloatType>(z.real()),
+                                   static_cast<ToFloatType>(z.imag()));
+}
 
- template <typename ToFloatType, typename FromFloatType>
-  std::vector<std::complex<ToFloatType> > ConvertComplexVector(std::vector<std::complex<FromFloatType> > x) {
-    std::vector<std::complex<ToFloatType> > new_x;
-    for (auto element : x) {
-      new_x.push_back(std::complex<ToFloatType>(static_cast<ToFloatType>(element.real()),
-                                                static_cast<ToFloatType>(element.imag()) ) );
-    }
-    return new_x;
+template <typename ToFloatType, typename FromFloatType>
+std::vector<std::complex<ToFloatType> > ConvertComplexVector(
+    std::vector<std::complex<FromFloatType> > x) {
+  std::vector<std::complex<ToFloatType> > new_x;
+  for (auto element : x) {
+    new_x.push_back(
+        std::complex<ToFloatType>(static_cast<ToFloatType>(element.real()),
+                                  static_cast<ToFloatType>(element.imag())));
   }
+  return new_x;
+}
 
-  template <typename ToFloatType, typename FromFloatType>
-  std::vector<std::vector<std::complex<ToFloatType> > > ConvertComplexVectorVector(std::vector<std::vector<std::complex<FromFloatType> > > x) {
-    std::vector<std::vector<std::complex<ToFloatType> > > new_x;
-    std::vector<std::complex<ToFloatType> >  new_y;
-    for (auto y : x) {
-      new_y.clear();
-      for (auto element : y) {
-        new_y.push_back(std::complex<ToFloatType>(static_cast<ToFloatType>(element.real()),
-                                                  static_cast<ToFloatType>(element.imag()) ) );
-      }
-      new_x.push_back(new_y);
+template <typename ToFloatType, typename FromFloatType>
+std::vector<std::vector<std::complex<ToFloatType> > >
+ConvertComplexVectorVector(
+    std::vector<std::vector<std::complex<FromFloatType> > > x) {
+  std::vector<std::vector<std::complex<ToFloatType> > > new_x;
+  std::vector<std::complex<ToFloatType> > new_y;
+  for (auto y : x) {
+    new_y.clear();
+    for (auto element : y) {
+      new_y.push_back(
+          std::complex<ToFloatType>(static_cast<ToFloatType>(element.real()),
+                                    static_cast<ToFloatType>(element.imag())));
     }
-    return new_x;
+    new_x.push_back(new_y);
   }
+  return new_x;
+}
 
 }  // end of namespace nmie
 #endif  // SRC_NMIE_PRECISION_H_

src/nmie-web-impl.hpp

@@ -0,0 +1,79 @@
+#ifndef SRC_NMIE_WEB_IMPL_HPP_
+#define SRC_NMIE_WEB_IMPL_HPP_
+//**********************************************************************************
+//    Copyright (C) 2009-2024  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2024  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite at least one of the following references:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+//        calculation of electromagnetic near-field for a multilayered
+//        sphere," Computer Physics Communications, vol. 214, May 2017,
+//        pp. 225-230.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+//    @brief  Wrapper class around nMie function for ease of use
+//
+//**********************************************************************************
+#include <algorithm>
+#include <array>
+#include <cstdio>
+#include <cstdlib>
+#include <stdexcept>
+#include <vector>
+
+#include <emscripten/bind.h>
+#include <emscripten/val.h>
+#include "nmie-web.hpp"
+
+namespace nmie {
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+
+// from https://toughengineer.github.io/demo/dsp/fft-perf/
+template <typename FloatType = double>
+emscripten::val toJSFloat64Array(const std::vector<double>& v) {
+  emscripten::val view{emscripten::typed_memory_view(
+      v.size(), v.data())};  // make a view of local object
+
+  auto result = emscripten::val::global("Float64Array")
+                    .new_(v.size());  // make a JS typed array
+  result.call<void>(
+      "set",
+      view);  // set typed array values "on the JS side" using the memory view
+
+  return result;
+}
+
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+template <typename FloatType>
+emscripten::val MultiLayerMieWeb<FloatType>::GetFieldEabs() {
+  auto Eabs = this->MultiLayerMie<FloatType>::GetFieldEabs();
+  return toJSFloat64Array(Eabs);
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+}  // end of namespace nmie
+#endif  // SRC_NMIE_WEB_IMPL_HPP_

src/nmie-web.hpp

@@ -0,0 +1,56 @@
+#ifndef SRC_NMIE_WEB_HPP_
+#define SRC_NMIE_WEB_HPP_
+//**********************************************************************************
+//    Copyright (C) 2009-2024  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2024  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite at least one of the following references:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+//        calculation of electromagnetic near-field for a multilayered
+//        sphere," Computer Physics Communications, vol. 214, May 2017,
+//        pp. 225-230.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//**********************************************************************************
+
+#include <array>
+#include <complex>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#include "nmie-applied-impl.hpp"
+
+#include <emscripten/bind.h>
+#include <emscripten/val.h>
+
+namespace nmie {
+
+template <typename FloatType = double>
+class MultiLayerMieWeb : public MultiLayerMieApplied<FloatType> {
+  // Will throw for any error!
+ public:
+  emscripten::val GetFieldEabs();
+
+};  // end of class MultiLayerMieWeb
+
+}  // end of namespace nmie
+#endif  // SRC_NMIE_WEB_HPP

tests/c++/go-speed-test.sh

@@ -15,7 +15,11 @@ rm -f $PROGRAM
 
 file=speed-test-applied.cc
 # g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
-g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+
+# with libtcmalloc
+# g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+
+g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -o $PROGRAM -march=native -mtune=native
 
 echo Should be:
 echo test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01

tests/c++/speed-test.cc

@@ -1,30 +1,34 @@
-//**********************************************************************************//
-//    Copyright (C) 2009-2016  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2016  Konstantin Ladutenko <kostyfisik@gmail.com>          //
-//                                                                                  //
-//    This file is part of scattnlay                                                //
-//                                                                                  //
-//    This program is free software: you can redistribute it and/or modify          //
-//    it under the terms of the GNU General Public License as published by          //
-//    the Free Software Foundation, either version 3 of the License, or             //
-//    (at your option) any later version.                                           //
-//                                                                                  //
-//    This program is distributed in the hope that it will be useful,               //
-//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
-//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
-//    GNU General Public License for more details.                                  //
-//                                                                                  //
-//    The only additional remark is that we expect that all publications            //
-//    describing work using this software, or all commercial products               //
-//    using it, cite the following reference:                                       //
-//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
-//        a multilayered sphere," Computer Physics Communications,                  //
-//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
-//                                                                                  //
-//    You should have received a copy of the GNU General Public License             //
-//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
-//**********************************************************************************//
+//**********************************************************************************
+//    Copyright (C) 2009-2016  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2016  Konstantin Ladutenko <kostyfisik@gmail.com>
+//
+//    This file is part of scattnlay
+//
+//    This program is free software: you can redistribute it and/or modify
+//    it under the terms of the GNU General Public License as published by
+//    the Free Software Foundation, either version 3 of the License, or
+//    (at your option) any later version.
+//
+//    This program is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//    GNU General Public License for more details.
+//
+//    The only additional remark is that we expect that all publications
+//    describing work using this software, or all commercial products
+//    using it, cite the following reference:
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+//        a multilayered sphere," Computer Physics Communications,
+//        vol. 180, Nov. 2009, pp. 2348-2354.
+//
+//    You should have received a copy of the GNU General Public License
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//**********************************************************************************
 
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
 #include <algorithm>
 #include <complex>
 #include <functional>
@@ -32,47 +36,54 @@
 #include <stdexcept>
 #include <string>
 #include <vector>
-#include <stdlib.h>
-#include <stdio.h>
-#include <time.h>
-#include <string.h>
-//sudo aptitude install libgoogle-perftools-dev
-//#include <google/heap-profiler.h>
-#include "../../src/nmie.hpp"
-//#include "../../src/nmie-precision.hpp"
-//#include "../../src/nmie-impl.hpp"
+// sudo aptitude install libgoogle-perftools-dev
+// #include <google/heap-profiler.h>
+#include "../../src/nmie-applied-impl.hpp"
 
 timespec diff(timespec start, timespec end);
-const double PI=3.14159265358979323846;
-template<class T> inline T pow2(const T value) {return value*value;}
-
+const double PI = 3.14159265358979323846;
+template <class T>
+inline T pow2(const T value) {
+  return value * value;
+}
 
-//***********************************************************************************//
-// This is the main function of 'scattnlay', here we read the parameters as          //
-// arguments passed to the program which should be executed with the following       //
-// syntaxis:                                                                         //
-// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]  //
-//                                                                                   //
-// When all the parameters were correctly passed we setup the integer L (the         //
-// number of layers) and the arrays x and m, containing the size parameters and      //
-// refractive indexes of the layers, respectively and call the function nMie.        //
-// If the calculation is successful the results are printed with the following       //
-// format:                                                                           //
-//                                                                                   //
-//    * If no comment was passed:                                                    //
-//        'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                                    //
-//                                                                                   //
-//    * If a comment was passed:                                                     //
-//        'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                           //
-//***********************************************************************************//
-int main(int argc, char *argv[]) {
+//***********************************************************************************
+// This is the main function of 'scattnlay', here we read the parameters as
+// arguments passed to the program which should be executed with the following
+// syntaxis:
+// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c
+// comment]
+//
+// When all the parameters were correctly passed we setup the integer L (the
+// number of layers) and the arrays x and m, containing the size parameters and
+// refractive indexes of the layers, respectively and call the function nMie.
+// If the calculation is successful the results are printed with the following
+// format:
+//
+//    * If no comment was passed:
+//        'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'
+//
+//    * If a comment was passed:
+//        'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'
+//***********************************************************************************
+int main(int argc, char* argv[]) {
   try {
     std::vector<std::string> args;
     args.assign(argv, argv + argc);
-    std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + args[0]
-			  + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] "
-			  + "[-t ti tf nt] [-c comment]\n");
-    enum mode_states {read_L, read_x, read_mr, read_mi, read_ti, read_tf, read_nt, read_comment};
+    std::string error_msg(std::string("Insufficient parameters.\nUsage: ") +
+                          args[0] +
+                          " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] " +
+                          "[-t ti tf nt] [-c comment]\n");
+    enum mode_states {
+      read_L,
+      read_x,
+      read_mr,
+      read_mi,
+      read_ti,
+      read_tf,
+      read_nt,
+      read_comment
+    };
     // for (auto arg : args) std::cout<< arg <<std::endl;
     std::string comment;
     int has_comment = 0;
@@ -86,10 +97,11 @@ int main(int argc, char *argv[]) {
 
     double ti = 0.0, tf = 90.0;
     int nt = 0;
-    if (argc < 5) throw std::invalid_argument(error_msg);
+    if (argc < 5)
+      throw std::invalid_argument(error_msg);
 
-    //strcpy(comment, "");
-    // for (i = 1; i < argc; i++) {
+    // strcpy(comment, "");
+    //  for (i = 1; i < argc; i++) {
     int mode = -1;
     double tmp_mr;
     for (auto arg : args) {
@@ -99,84 +111,85 @@ int main(int argc, char *argv[]) {
 
       // Detecting new read mode (if it is a valid -key)
       if (arg == "-l") {
-	mode = read_L;
-	continue;
+        mode = read_L;
+        continue;
       }
       if (arg == "-t") {
-	if ((mode != read_x) && (mode != read_comment))
-	  throw std::invalid_argument(std::string("Unfinished layer!\n")
-							 +error_msg);
-	mode = read_ti;
-	continue;
+        if ((mode != read_x) && (mode != read_comment))
+          throw std::invalid_argument(std::string("Unfinished layer!\n") +
+                                      error_msg);
+        mode = read_ti;
+        continue;
       }
       if (arg == "-c") {
-	if ((mode != read_x) && (mode != read_nt))
-	  throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg);
-	mode = read_comment;
-	continue;
+        if ((mode != read_x) && (mode != read_nt))
+          throw std::invalid_argument(
+              std::string("Unfinished layer or theta!\n") + error_msg);
+        mode = read_comment;
+        continue;
       }
       // Reading data. For invalid date the exception will be thrown
       // with the std:: and catched in the end.
       if (mode == read_L) {
-	L = std::stoi(arg);
-	mode = read_x;
-	continue;
+        L = std::stoi(arg);
+        mode = read_x;
+        continue;
       }
       if (mode == read_x) {
-	x.push_back(std::stod(arg));
-	mode = read_mr;
-	continue;
+        x.push_back(std::stod(arg));
+        mode = read_mr;
+        continue;
       }
       if (mode == read_mr) {
-	tmp_mr = std::stod(arg);
-	mode = read_mi;
-	continue;
+        tmp_mr = std::stod(arg);
+        mode = read_mi;
+        continue;
       }
       if (mode == read_mi) {
-	m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
-	mode = read_x;
-	continue;
+        m.push_back(std::complex<double>(tmp_mr, std::stod(arg)));
+        mode = read_x;
+        continue;
       }
       if (mode == read_ti) {
-	ti = std::stod(arg);
-	mode = read_tf;
-	continue;
+        ti = std::stod(arg);
+        mode = read_tf;
+        continue;
       }
       if (mode == read_tf) {
-	tf = std::stod(arg);
-	mode = read_nt;
-	continue;
+        tf = std::stod(arg);
+        mode = read_nt;
+        continue;
       }
       if (mode == read_nt) {
-	nt = std::stoi(arg);
+        nt = std::stoi(arg);
         Theta.resize(nt);
         S1.resize(nt);
         S2.resize(nt);
         S1w.resize(nt);
         S2w.resize(nt);
-	continue;
+        continue;
       }
-      if (mode ==  read_comment) {
-	comment = arg;
+      if (mode == read_comment) {
+        comment = arg;
         has_comment = 1;
-	continue;
+        continue;
       }
     }
-    if ( (x.size() != m.size()) || (L != x.size()) )
-      throw std::invalid_argument(std::string("Broken structure!\n")
-							 +error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) )
-      throw std::invalid_argument(std::string("Empty structure!\n")
-							 +error_msg);
+    if ((x.size() != m.size()) || (L != x.size()))
+      throw std::invalid_argument(std::string("Broken structure!\n") +
+                                  error_msg);
+    if ((0 == m.size()) || (0 == x.size()))
+      throw std::invalid_argument(std::string("Empty structure!\n") +
+                                  error_msg);
 
     if (nt < 0) {
       printf("Error reading Theta.\n");
       return -1;
     } else if (nt == 1) {
-      Theta[0] = ti*PI/180.0;
+      Theta[0] = ti * PI / 180.0;
     } else {
       for (i = 0; i < nt; i++) {
-      Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
+        Theta[i] = (ti + (double)i * (tf - ti) / (nt - 1)) * PI / 180.0;
       }
     }
     timespec time1, time2;
@@ -184,59 +197,60 @@ int main(int argc, char *argv[]) {
     long ctime_nsec, best_c;
     long cpptime_sec, ctime_sec;
     long repeats = 150;
-    //HeapProfilerStart("heapprof");
+    // HeapProfilerStart("heapprof");
     printf("Start\n");
     do {
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
-      for (int i = 0; i<repeats; ++i) {
-    	nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw,
-    			   &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
-	// break;
+      for (int i = 0; i < repeats; ++i) {
+        nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw,
+                   &gw, &Albedow, S1w, S2w);
+        // break;
       }
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
-      cpptime_nsec = diff(time1,time2).tv_nsec;
-      cpptime_sec = diff(time1,time2).tv_sec;
+      cpptime_nsec = diff(time1, time2).tv_nsec;
+      cpptime_sec = diff(time1, time2).tv_sec;
 
-      printf("-- C++ time consumed %lg sec\n", cpptime_sec+(cpptime_nsec/1e9));
+      printf("-- C++ time consumed %lg sec\n",
+             cpptime_sec + (cpptime_nsec / 1e9));
       repeats *= 10;
       // break;
     } while (cpptime_sec < 3 && ctime_sec < 3);
 
-        printf("Finish\n");
+    printf("Finish\n");
 
     // if (has_comment) {
-    //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n",
+    //   comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     // } else {
-    //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw,
+    //   Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     // }
 
     // if (nt > 0) {
-    //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+    //   printf(" Theta,         S1.r,         S1.i,         S2.r, S2.i\n");
 
     //   for (i = 0; i < nt; i++) {
-    //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
+    //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI,
+    //     S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
     //   }
     // }
 
-  } catch( const std::invalid_argument &ia ) {
+  } catch (const std::invalid_argument& ia) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
   }
-    return 0;
+  return 0;
 }
 
-
-
-timespec diff(timespec start, timespec end)
-{
-	timespec temp;
-	if ((end.tv_nsec-start.tv_nsec)<0) {
-		temp.tv_sec = end.tv_sec-start.tv_sec-1;
-		temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
-	} else {
-		temp.tv_sec = end.tv_sec-start.tv_sec;
-		temp.tv_nsec = end.tv_nsec-start.tv_nsec;
-	}
-	return temp;
+timespec diff(timespec start, timespec end) {
+  timespec temp;
+  if ((end.tv_nsec - start.tv_nsec) < 0) {
+    temp.tv_sec = end.tv_sec - start.tv_sec - 1;
+    temp.tv_nsec = 1000000000 + end.tv_nsec - start.tv_nsec;
+  } else {
+    temp.tv_sec = end.tv_sec - start.tv_sec;
+    temp.tv_nsec = end.tv_nsec - start.tv_nsec;
+  }
+  return temp;
 }

tests/test_py.py

@@ -56,14 +56,15 @@ class TestStringMethods(unittest.TestCase):
                 continue
             print("Using solver: ", solver)
             for case in test_cases:
-                print("test case:", case)
                 solver.SetLayersSize(case[0])
                 solver.SetLayersIndex(case[1])
                 solver.RunMieCalculation()
                 Qext = solver.GetQext()
                 Qsca = solver.GetQsca()
-                print("ext tol:", (case[2] - Qext) / Qext)
-                print("sca tol:", (case[3] - Qsca) / Qsca)
+                if case == test_cases[0]:
+                    print("test case:", case)
+                    print("ext tol:", (case[2] - Qext) / Qext)
+                    print("sca tol:", (case[3] - Qsca) / Qsca)
                 self.assertTrue((case[2] - Qext) / Qext < tol)
                 self.assertTrue((case[3] - Qsca) / Qsca < tol)