calc_ab for mesomie compiles

src/mesomie.hpp

@@ -1,5 +1,5 @@
-#ifndef SRC_NMIE_BASIC_HPP_
-#define SRC_NMIE_BASIC_HPP_
+#ifndef SRC_MESOMIE_HPP_
+#define SRC_MESOMIE_HPP_
 //******************************************************************************
 //    Copyright (C) 2009-2022  Ovidio Pena <ovidio@bytesfall.com>
 //    Copyright (C) 2013-2022  Konstantin Ladutenko <kostyfisik@gmail.com>
@@ -48,10 +48,6 @@
 //
 // Hereinafter all equations numbers refer to [2]
 //******************************************************************************
-/*******************************************************************************
- * sdfa
- * sdf
- */
 #include <iomanip>
 #include <iostream>
 #include <stdexcept>
@@ -73,7 +69,6 @@ void MesoMie<FloatType>::calc_ab(int nmax,
                                  FloatType d_perp) {
   an_.resize(nmax, static_cast<FloatType>(0.0));
   bn_.resize(nmax, static_cast<FloatType>(0.0));
-
   std::vector<std::complex<FloatType>>      //
       D1_xd(nmax + 1), D3_xd(nmax + 1),     //
       D1_xm(nmax + 1), D3_xm(nmax + 1),     //
@@ -82,6 +77,39 @@ void MesoMie<FloatType>::calc_ab(int nmax,
 
   evalPsiZetaD1D3(std::complex<FloatType>(xd), Psi_xd, Zeta_xd, D1_xd, D3_xd);
   evalPsiZetaD1D3(std::complex<FloatType>(xm), Psi_xm, Zeta_xm, D1_xm, D3_xm);
+
+  for (int n = 1; n <= nmax; n++) {
+    an_[n] = Psi_xd[n] *
+             (                                                       //
+                 eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] +     //
+                 (eps_m - eps_d) *                                   //
+                     (                                               //
+                         n * (n + 1) * d_perp +                      //
+                         xd * D1_xd[n] * xm * D1_xm[n] * d_parallel  //
+                         ) /                                         //
+                     R                                               //
+                 ) /                                                 //
+             Zeta_xd[n] *
+             (                                                       //
+                 eps_m * xd * D3_xd[n] - eps_d * xm * D1_xm[n] +     //
+                 (eps_m - eps_d) *                                   //
+                     (                                               //
+                         n * (n + 1) * d_perp +                      //
+                         xd * D3_xd[n] * xm * D1_xm[n] * d_parallel  //
+                         ) /                                         //
+                     R                                               //
+             );
+    bn_[n] = Psi_xd[n] *
+             (                                         //
+                 xd * D1_xd[n] - xm * D1_xm[n] +       //
+                 (xm * xm - xd * xd) * d_parallel / R  //
+                 ) /                                   //
+             Zeta_xd[n] *
+             (                                         //
+                 xd * D3_xd[n] - xm * D1_xm[n] +       //
+                 (xm * xm - xd * xd) * d_parallel / R  //
+             );                                        //
+  }
 }
 
 }  // namespace nmie

src/nearfield.cc

@@ -1,69 +1,84 @@
-//**********************************************************************************//
-//    Copyright (C) 2009-2022  Ovidio Pena <ovidio@bytesfall.com>                   //
-//    Copyright (C) 2013-2022  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-2022  Ovidio Pena <ovidio@bytesfall.com>
+//    Copyright (C) 2013-2022  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 <complex>
+#include <cstdio>
 #include <iostream>
 #include <stdexcept>
 #include <string>
 #include <vector>
-#include <cstdio>
 
 #include "nmie.hpp"
 
-const double PI=3.14159265358979323846;
-
-//***********************************************************************************//
-// 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:                                                                         //
-// ./fieldnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...]                             //
-//             -p xi xf nx yi yf ny zi zf nz [-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:                                                                           //
-//                                                                                   //
-// 'X, Y, Z, Ex.r, Ex.i, Ey.r, Ey.i, Ez.r, Ez.i, Hx.r, Hx.i, Hy.r, Hy.i, Hz.r, Hz.i' //
-//                                                                                   //
-//***********************************************************************************//
-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:
+// ./fieldnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...]
+//             -p xi xf nx yi yf ny zi zf nz [-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:
+//
+// 'X, Y, Z, Ex.r, Ex.i, Ey.r, Ey.i, Ez.r, Ez.i, Hx.r, Hx.i, Hy.r, Hy.i, Hz.r,
+// Hz.i'
+//
+//******************************************************************************
+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 ...] "
-                          + " -p xi xf nx yi yf ny zi zf nz [-c comment]\n");
-    enum mode_states {read_L, read_x, read_mr, read_mi, read_xi, read_xf, read_nx, read_yi, read_yf, read_ny, read_zi, read_zf, read_nz, 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 ...] " +
+                          " -p xi xf nx yi yf ny zi zf nz [-c comment]\n");
+    enum mode_states {
+      read_L,
+      read_x,
+      read_mr,
+      read_mi,
+      read_xi,
+      read_xf,
+      read_nx,
+      read_yi,
+      read_yf,
+      read_ny,
+      read_zi,
+      read_zf,
+      read_nz,
+      read_comment
+    };
     // for (auto arg : args) std::cout<< arg <<std::endl;
     std::string comment;
     int has_comment = 0;
@@ -76,11 +91,12 @@ int main(int argc, char *argv[]) {
     double dx = 0.0, dy = 0.0, dz = 0.0;
     int nx = 0, ny = 0, nz = 0;
     long total_points = 0;
-    if (argc < 5) throw std::invalid_argument(error_msg);
+    if (argc < 5)
+      throw std::invalid_argument(error_msg);
 
     int mode = -1;
     double tmp_mr;
-    for (const auto &arg : args) {
+    for (const auto& arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
@@ -93,14 +109,16 @@ int main(int argc, char *argv[]) {
 
       if (arg == "-p") {
         if ((mode != read_x) && (mode != read_comment))
-          throw std::invalid_argument(std::string("Unfinished layer!\n") + error_msg);
+          throw std::invalid_argument(std::string("Unfinished layer!\n") +
+                                      error_msg);
         mode = read_xi;
         continue;
       }
 
       if (arg == "-c") {
         if ((mode != read_x) && (mode != read_nz))
-          throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg);
+          throw std::invalid_argument(
+              std::string("Unfinished layer or theta!\n") + error_msg);
         mode = read_comment;
         continue;
       }
@@ -126,7 +144,7 @@ int main(int argc, char *argv[]) {
       }
 
       if (mode == read_mi) {
-        m.emplace_back( tmp_mr,std::stod(arg) );
+        m.emplace_back(tmp_mr, std::stod(arg));
         mode = read_x;
         continue;
       }
@@ -181,9 +199,12 @@ int main(int argc, char *argv[]) {
 
       if (mode == read_nz) {
         nz = std::stoi(arg);
-        total_points = nx*ny*nz;
+        total_points = nx * ny * nz;
         if (total_points <= 0)
-          throw std::invalid_argument(std::string("Nothing to do! You must define the grid to calculate the fields.\n") + error_msg);
+          throw std::invalid_argument(
+              std::string("Nothing to do! You must define the grid to "
+                          "calculate the fields.\n") +
+              error_msg);
 
         Xp.resize(total_points);
         Yp.resize(total_points);
@@ -198,66 +219,73 @@ int main(int argc, char *argv[]) {
         continue;
       }
 
-      if (mode ==  read_comment) {
+      if (mode == read_comment) {
         comment = arg;
         has_comment = 1;
         continue;
       }
     }
 
-    if ( (x.size() != m.size()) || (L != x.size()) )
-      throw std::invalid_argument(std::string("Broken structure!\n") + error_msg);
-    if ( (m.empty()) || ( x.empty()) )
-      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 ((m.empty()) || (x.empty()))
+      throw std::invalid_argument(std::string("Empty structure!\n") +
+                                  error_msg);
 
     if (nx == 1)
       dx = 0.0;
     else
-      dx = (xf - xi)/(nx - 1);
+      dx = (xf - xi) / (nx - 1);
 
     if (ny == 1)
       dy = 0.0;
     else
-      dy = (yf - yi)/(ny - 1);
+      dy = (yf - yi) / (ny - 1);
 
     if (nz == 1)
       dz = 0.0;
     else
-      dz = (zf - zi)/(nz - 1);
+      dz = (zf - zi) / (nz - 1);
 
     for (int i = 0; i < nx; i++) {
       for (int j = 0; j < ny; j++) {
         for (int k = 0; k < nz; k++) {
-          Xp[i*ny*nz + j*nz + k] = xi + (double)i*dx;
-          Yp[i*ny*nz + j*nz + k] = yi + (double)j*dy;
-          Zp[i*ny*nz + j*nz + k] = zi + (double)k*dz;
+          Xp[i * ny * nz + j * nz + k] = xi + (double)i * dx;
+          Yp[i * ny * nz + j * nz + k] = yi + (double)j * dy;
+          Zp[i * ny * nz + j * nz + k] = zi + (double)k * dz;
         }
       }
     }
 
-    int nmax = nmie::nField(L, -1, x, m, -1, nmie::Modes::kAll, nmie::Modes::kAll, total_points, Xp, Yp, Zp, E, H);
+    int nmax = nmie::nField(L, -1, x, m, -1, nmie::Modes::kAll,
+                            nmie::Modes::kAll, total_points, Xp, Yp, Zp, E, H);
     printf("Number of multipoles used in Mie series nmax=%i\n", nmax);
     if (has_comment)
       printf("%6s\n", comment.c_str());
 
     if (total_points > 0) {
-      printf("         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i\n");
+      printf(
+          "         X,          Y,          Z,         Ex.r,         Ex.i,     "
+          "    Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,   "
+          "      Hx.i,         Hy.r,         Hy.i,         Hz.r,         "
+          "Hz.i\n");
 
       for (long i = 0; i < total_points; i++) {
-        printf("%10.7f, %10.7f, %10.7f, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n",
-               Xp[i], Yp[i], Zp[i],
-               E[i][0].real(), E[i][0].imag(), E[i][1].real(), E[i][1].imag(), E[i][2].real(), E[i][2].imag(),
-               H[i][0].real(), H[i][0].imag(), H[i][1].real(), H[i][1].imag(), H[i][2].real(), H[i][2].imag());
+        printf(
+            "%10.7f, %10.7f, %10.7f, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, "
+            "%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n",
+            Xp[i], Yp[i], Zp[i], E[i][0].real(), E[i][0].imag(), E[i][1].real(),
+            E[i][1].imag(), E[i][2].real(), E[i][2].imag(), H[i][0].real(),
+            H[i][0].imag(), H[i][1].real(), H[i][1].imag(), H[i][2].real(),
+            H[i][2].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;
 }
-
-

tests/test_bulk_sphere.cc

@@ -1,80 +1,143 @@
-#include "gtest/gtest.h"
 #include "../src/nmie-basic.hpp"
+#include "../src/mesomie.hpp"
+#include "gtest/gtest.h"
 
-// TODO fails for MP with 100 digits. And 16 digits, which should be equal to double precision.
+// TODO fails for MP with 100 digits. And 16 digits, which should be equal to
+// double precision.
 #ifndef MULTI_PRECISION
-//TEST(BulkSphere, DISABLED_ArgPi) {
+// TEST(BulkSphere, DISABLED_ArgPi) {
+//******************************************************************************
 TEST(BulkSphere, ArgPi) {
-  std::vector<double> WLs{50, 80, 100,200, 400}; //nm
+  std::vector<double> WLs{50, 80, 100, 200, 400};  // nm
   double host_index = 2.;
-  double core_radius = 100.; //nm
+  double core_radius = 100.;  // nm
   double delta = 1e-5;
   nmie::MultiLayerMie<nmie::FloatType> nmie;
-  nmie.SetLayersIndex({std::complex<double>(4,0)});
-  for (auto WL:WLs) {
-    nmie.SetLayersSize({2*nmie::PI_*host_index*core_radius/(WL+delta)});
+  nmie.SetLayersIndex({std::complex<double>(4, 0)});
+  for (auto WL : WLs) {
+    nmie.SetLayersSize(
+        {2 * nmie::PI_ * host_index * core_radius / (WL + delta)});
     nmie.RunMieCalculation();
     double Qabs_p = std::abs(static_cast<double>(nmie.GetQabs()));
 
-    nmie.SetLayersSize({2*nmie::PI_*host_index*core_radius/(WL-delta)});
+    nmie.SetLayersSize(
+        {2 * nmie::PI_ * host_index * core_radius / (WL - delta)});
     nmie.RunMieCalculation();
     double Qabs_m = std::abs(static_cast<double>(nmie.GetQabs()));
 
-    nmie.SetLayersSize({2*nmie::PI_*host_index*core_radius/(WL)});
+    nmie.SetLayersSize({2 * nmie::PI_ * host_index * core_radius / (WL)});
     nmie.RunMieCalculation();
     double Qabs = std::abs(static_cast<double>(nmie.GetQabs()));
-    EXPECT_GT(Qabs_p+Qabs_m, Qabs);
+    EXPECT_GT(Qabs_p + Qabs_m, Qabs);
   }
 }
 #endif
 
-//TEST(BulkSphere, DISABLED_HandlesInput) {
+//******************************************************************************
+// TEST(BulkSphere, DISABLED_HandlesInput) {
 TEST(BulkSphere, HandlesInput) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
   // A list of tests for a bulk sphere from
   // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
   // table 1: sphere size and refractive index
   // followed by resulting extinction and scattering efficiencies
-  std::vector< std::tuple< double, std::complex<double>, double, double, char > >
-      parameters_and_results
-      {
+  std::vector<std::tuple<double, std::complex<double>, double, double, char> >
+      parameters_and_results{
           // x, {Re(m), Im(m)}, Qext, Qsca, test_name
-          {0.099, {0.75,0}, 7.417859e-06, 7.417859e-06, 'a'},
-          {0.101, {0.75,0}, 8.033538e-06, 8.033538e-06, 'b'},
-          {10,    {0.75,0},     2.232265, 2.232265, 'c'},
-          {100,   {1.33,1e-5}, 2.101321, 2.096594, 'e'},
-          {0.055, {1.5, 1},    0.10149104, 1.131687e-05, 'g'},
-          {0.056, {1.5, 1},   0.1033467, 1.216311e-05, 'h'},
-          {100,   {1.5, 1},    2.097502, 1.283697, 'i'},
-          {1,     {10,  10},   2.532993, 2.049405, 'k'},
-          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
-          {100,   {10,  10,},  2.071124, 1.836785, 'l'},
-          {10000, {1.33,1e-5}, 2.004089, 1.723857, 'f'},
-          {10000, {1.5, 1},    2.004368, 1.236574, 'j'},
-          {10000, {10,  10},   2.005914, 1.795393, 'm'},
+          {0.099, {0.75, 0}, 7.417859e-06, 7.417859e-06, 'a'},
+          {0.101, {0.75, 0}, 8.033538e-06, 8.033538e-06, 'b'},
+          {10, {0.75, 0}, 2.232265, 2.232265, 'c'},
+          {100, {1.33, 1e-5}, 2.101321, 2.096594, 'e'},
+          {0.055, {1.5, 1}, 0.10149104, 1.131687e-05, 'g'},
+          {0.056, {1.5, 1}, 0.1033467, 1.216311e-05, 'h'},
+          {100, {1.5, 1}, 2.097502, 1.283697, 'i'},
+          {1, {10, 10}, 2.532993, 2.049405, 'k'},
+          {1000, {0.75, 0}, 1.997908, 1.997908, 'd'},
+          {100,
+           {
+               10,
+               10,
+           },
+           2.071124,
+           1.836785,
+           'l'},
+          {10000, {1.33, 1e-5}, 2.004089, 1.723857, 'f'},
+          {10000, {1.5, 1}, 2.004368, 1.236574, 'j'},
+          {10000, {10, 10}, 2.005914, 1.795393, 'm'},
       };
-  for (const auto &data : parameters_and_results) {
+  for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
     auto m = std::get<1>(data);
-//    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
+    //    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
 
     nmie.SetLayersSize({x});
     nmie.SetLayersIndex({m});
-//    nmie.SetMaxTerms(Nstop);
+    //    nmie.SetMaxTerms(Nstop);
     nmie.RunMieCalculation();
     double Qext = static_cast<double>(nmie.GetQext());
     double Qsca = static_cast<double>(nmie.GetQsca());
-    double Qext_Du =  std::get<2>(data);
-    double Qsca_Du =  std::get<3>(data);
+    double Qext_Du = std::get<2>(data);
+    double Qsca_Du = std::get<3>(data);
     EXPECT_FLOAT_EQ(Qext_Du, Qext)
-              << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
-              << "\nnmax_ = " << nmie.GetMaxTerms();
+        << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
+        << "\nnmax_ = " << nmie.GetMaxTerms();
     EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
-              << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
+        << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
   }
+}
 
+//******************************************************************************
+TEST(BulkSphere, MesoMie) {
+  nmie::MultiLayerMie<nmie::FloatType> nmie;
+  // A list of tests for a bulk sphere from
+  // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
+  // table 1: sphere size and refractive index
+  // followed by resulting extinction and scattering efficiencies
+  std::vector<std::tuple<double, std::complex<double>, double, double, char> >
+      parameters_and_results{
+          // x, {Re(m), Im(m)}, Qext, Qsca, test_name
+          {0.099, {0.75, 0}, 7.417859e-06, 7.417859e-06, 'a'},
+          {0.101, {0.75, 0}, 8.033538e-06, 8.033538e-06, 'b'},
+          {10, {0.75, 0}, 2.232265, 2.232265, 'c'},
+          {100, {1.33, 1e-5}, 2.101321, 2.096594, 'e'},
+          {0.055, {1.5, 1}, 0.10149104, 1.131687e-05, 'g'},
+          {0.056, {1.5, 1}, 0.1033467, 1.216311e-05, 'h'},
+          {100, {1.5, 1}, 2.097502, 1.283697, 'i'},
+          {1, {10, 10}, 2.532993, 2.049405, 'k'},
+          {1000, {0.75, 0}, 1.997908, 1.997908, 'd'},
+          {100,
+           {
+               10,
+               10,
+           },
+           2.071124,
+           1.836785,
+           'l'},
+          {10000, {1.33, 1e-5}, 2.004089, 1.723857, 'f'},
+          {10000, {1.5, 1}, 2.004368, 1.236574, 'j'},
+          {10000, {10, 10}, 2.005914, 1.795393, 'm'},
+      };
+  for (const auto& data : parameters_and_results) {
+    auto x = std::get<0>(data);
+    auto m = std::get<1>(data);
+    //    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
+
+    nmie.SetLayersSize({x});
+    nmie.SetLayersIndex({m});
+    //    nmie.SetMaxTerms(Nstop);
+    nmie.RunMieCalculation();
+    double Qext = static_cast<double>(nmie.GetQext());
+    double Qsca = static_cast<double>(nmie.GetQsca());
+    double Qext_Du = std::get<2>(data);
+    double Qsca_Du = std::get<3>(data);
+    EXPECT_FLOAT_EQ(Qext_Du, Qext)
+        << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
+        << "\nnmax_ = " << nmie.GetMaxTerms();
+    EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
+        << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
+  }
 }
-int main(int argc, char **argv) {
+int main(int argc, char** argv) {
   testing::InitGoogleTest(&argc, argv);
   return RUN_ALL_TESTS();
 }