Preparing scattnlay to use pybind11 instead of Cython for Python bindings (not yet working).

README.md

@@ -1,6 +1,6 @@
 ![output example](/doc/OutputExample.png)
 **Output example:** Field distribution inside layered Si\Ag\Si sphere
-and Poynting vector distribution in Ag sphere with poweflow lines
+and Poynting vector distribution in Ag sphere with powerflow lines
 calculated with Scattnlay (scripts  field-SiAgSi-flow.py and
 field-Ag-flow.py from example section as [revision](https://github.com/ovidiopr/scattnlay/commit/57c7261705a5776f78420c1f486e929517d5f584) ).
 

setup_pb11.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2018 Ovidio Peña Rodríguez <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. Peña 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. Peña-Rodríguez, "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/>.
+
+__version__ = '2.2'
+__title__ = 'Calculation of the scattering of EM radiation by a multilayered sphere'
+__mod__ = 'python-scattnlay'
+__author__ = 'Ovidio Peña Rodríguez'
+__email__ = 'ovidio@bytesfall.com'
+__url__ = 'https://github.com/ovidiopr/scattnlay'
+__download_url__ = 'https://github.com/ovidiopr/scattnlay/archive/v2.2.0.tar.gz'
+
+from distutils.core import setup
+from distutils.extension import Extension
+import numpy as np
+import pybind11 as pb
+
+setup(name = __mod__,
+      version = __version__,
+      description = __title__,
+      long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
+radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
+resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
+considerable improvement over previous implementations of similar algorithms. For details see: \
+O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
+      author = __author__,
+      author_email = __email__,
+      maintainer = __author__,
+      maintainer_email = __email__,
+      keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
+      url = __url__,
+      download_url = __download_url__,
+      license = 'GPL',
+      platforms = 'any',
+      ext_modules = [Extension("scattnlay",
+                               ["src/nmie.cc", "src/pb11_nmie.cc", "src/pb11_wrapper.cc"],
+                               language = "c++",
+                               include_dirs = [np.get_include(), pb.get_include(True)], 
+                               extra_compile_args=['-std=c++11']),
+                     Extension("scattnlay_mp",
+                               ["src/nmie.cc", "src/pb11_nmie.cc", "src/pb11_wrapper_mp.cc"],
+                               language = "c++",
+                               include_dirs = [np.get_include(), pb.get_include(True)], 
+                               extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
+)
+

src/nmie-impl.hpp

@@ -94,9 +94,7 @@ namespace nmie {
     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())
-						)
-		      );
+                                                static_cast<ToFloatType>(element.imag()) ) );
     }
     return new_x;
   }
@@ -108,10 +106,8 @@ namespace nmie {
     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_y.push_back(std::complex<ToFloatType>(static_cast<ToFloatType>(element.real()),
+                                                  static_cast<ToFloatType>(element.imag()) ) );
       }
       new_x.push_back(new_y);
     }
@@ -391,8 +387,8 @@ namespace nmie {
     std::complex<FloatType> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
     std::complex<FloatType> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
     // std::cout<< std::setprecision(100)
-    // 	     << "Ql "	<< PsiXL
-    // 	     <<std::endl;
+    //          << "Ql "        << PsiXL
+    //          << std::endl;
 
 
     return Num/Denom;
@@ -467,7 +463,7 @@ namespace nmie {
 
     // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
     PsiZeta_[0] = static_cast<FloatType>(0.5)*(static_cast<FloatType>(1.0) - std::complex<FloatType>(nmm::cos(2.0*z.real()), nmm::sin(2.0*z.real()))
-		       *static_cast<FloatType>(nmm::exp(-2.0*z.imag())));
+                 *static_cast<FloatType>(nmm::exp(-2.0*z.imag())));
     D3[0] = std::complex<FloatType>(0.0, 1.0);
 
     for (int n = 1; n <= nmax_; n++) {
@@ -701,7 +697,7 @@ namespace nmie {
       //*************************************************//
       // Upward recurrence for Q - equations (19a) and (19b)
       Num = std::complex<FloatType>(nmm::exp(-2.0*(z1.imag() - z2.imag())), 0.0)
-	*std::complex<FloatType>(nmm::cos(-2.0*z2.real()) - nmm::exp(-2.0*z2.imag()), nmm::sin(-2.0*z2.real()));
+      *std::complex<FloatType>(nmm::cos(-2.0*z2.real()) - nmm::exp(-2.0*z2.imag()), nmm::sin(-2.0*z2.real()));
       Denom = std::complex<FloatType>(nmm::cos(-2.0*z1.real()) - nmm::exp(-2.0*z1.imag()), nmm::sin(-2.0*z1.real()));
       Q[l][0] = Num/Denom;
 
@@ -984,17 +980,17 @@ namespace nmie {
       if (cabs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
       else {
         //throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
-	std::cout<< std::setprecision(100)
-		 << "Warning: Potentially unstable calculation of aln[0]["
-		 << n << "] = "<< aln_[0][n] <<std::endl;
+        std::cout<< std::setprecision(100)
+                 << "Warning: Potentially unstable calculation of aln[0]["
+                 << n << "] = "<< aln_[0][n] <<std::endl;
         aln_[0][n] = 0.0;
       }
       if (cabs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
       else {
         //throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
-	std::cout<< std::setprecision(100)
-		 << "Warning: Potentially unstable calculation of bln[0]["
-		 << n << "] = "<< bln_[0][n] <<std::endl;
+        std::cout<< std::setprecision(100)
+                 << "Warning: Potentially unstable calculation of bln[0]["
+                 << n << "] = "<< bln_[0][n] <<std::endl;
         bln_[0][n] = 0.0;
       }
     }
@@ -1065,7 +1061,7 @@ namespace nmie {
 
       // Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
       std::complex<FloatType> En = ipow[n1 % 4]
-	*static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
+      *static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
       for (int i = 0; i < 3; i++) {
         // electric field E [V m - 1] = EF*E0
         E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]

src/nmie.cc

@@ -78,7 +78,8 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
+  int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
+                  const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -133,7 +134,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(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 nMie(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!");
@@ -203,7 +207,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(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 nMie(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::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -235,7 +242,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(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 nMie(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::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -269,7 +279,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(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) {
+  int nMie(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::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -296,7 +309,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nField(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const unsigned int ncoord, const std::vector<double>& Xp_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
+  int nField(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m,
+             const int nmax, const unsigned int ncoord,
+             const std::vector<double>& Xp_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec,
+             std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
     if (x.size() != L || m.size() != L)
       throw std::invalid_argument("Declared number of layers do not fit x and m!");
     if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord

src/pb11_nmie.cc

@@ -0,0 +1,174 @@
+//**********************************************************************************//
+//    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 <complex>
+#include <vector>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include "nmie.hpp"
+//#include "pb11_nmie.hpp"
+
+
+namespace py = pybind11;
+
+
+py::tuple scattcoeffs(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                      py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                      int nmax, int pl)
+{
+  if (x.ndim() != 2)
+    throw std::runtime_error("The size parameter (x) should be 2-D NumPy array.");
+  if (m.ndim() != 2)
+    throw std::runtime_error("The relative refractive index (m) should be 2-D NumPy array.");
+
+
+  // allocate std::vector (to pass to the C++ function)
+  std::vector<std::vector<double> > x_cpp(x.size());
+  std::vector<std::vector<std::complex<double> > > m_cpp(m.size());
+
+  // copy py::array -> std::vector
+  std::memcpy(x_cpp.data(), x.data(), x.size()*sizeof(double));
+  std::memcpy(m_cpp.data(), m.data(), m.size()*sizeof(std::complex<double>));
+
+  // create std::vector (to get return from the C++ function)
+  std::vector<int> terms(x.shape(0));
+  std::vector<std::vector<std::complex<double> > > an, bn;
+  an.resize(x.shape(0));
+  bn.resize(x.shape(0));
+  
+  for (ssize_t i = 0; i < x.shape(0); i++) {
+    terms[i] = nmie::ScattCoeffs(x.shape(1), pl, x_cpp[i], m_cpp[i], nmax, an[i], bn[i]);
+  }
+
+  return py::make_tuple(terms, an, bn);
+
+}
+
+py::tuple scattnlay(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                    py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                    py::array_t<double, py::array::c_style | py::array::forcecast> theta,
+                    int nmax, int pl)
+{
+  if (x.ndim() != 2)
+    throw std::runtime_error("The size parameter (x) should be 2-D NumPy array.");
+  if (m.ndim() != 2)
+    throw std::runtime_error("The relative refractive index (m) should be 2-D NumPy array.");
+  if (theta.ndim() != 1)
+    throw std::runtime_error("The scattering angles (theta) should be 1-D NumPy array.");
+
+
+  // allocate std::vector (to pass to the C++ function)
+  std::vector<std::vector<double> > x_cpp(x.size());
+  std::vector<std::vector<std::complex<double> > > m_cpp(m.size());
+  std::vector<double> theta_cpp(theta.size());
+
+  // copy py::array -> std::vector
+  std::memcpy(x_cpp.data(), x.data(), x.size()*sizeof(double));
+  std::memcpy(m_cpp.data(), m.data(), m.size()*sizeof(std::complex<double>));
+  std::memcpy(theta_cpp.data(), theta.data(), theta.size()*sizeof(double));
+
+
+  // create std::vector (to get return from the C++ function)
+  std::vector<int> terms(x.shape(0));
+  std::vector<std::vector<std::complex<double> > > S1, S2;
+  S1.resize(x.shape(0));
+  S2.resize(x.shape(0));
+
+  std::vector<double> Qext(x.shape(0));
+  std::vector<double> Qsca(x.shape(0));
+  std::vector<double> Qabs(x.shape(0));
+  std::vector<double> Qbk(x.shape(0));
+  std::vector<double> Qpr(x.shape(0));
+  std::vector<double> g(x.shape(0));
+  std::vector<double> Albedo(x.shape(0));
+  
+  for (ssize_t i = 0; i < x.shape(0); i++) {
+    S1[i].resize(theta.shape(0));
+    S2[i].resize(theta.shape(0));
+
+    terms[i] = nmie::nMie(x.shape(1), pl, x_cpp[i], m_cpp[i], theta.shape(0), theta_cpp, nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], S1[i], S2[i]);
+  }
+
+  return py::make_tuple(terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+}
+
+py::tuple fieldnlay(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                    py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                    py::array_t<double, py::array::c_style | py::array::forcecast> coords,
+                    int nmax, int pl)
+{
+  if (x.ndim() != 2)
+    throw std::runtime_error("The size parameter (x) should be 2-D NumPy array.");
+  if (m.ndim() != 2)
+    throw std::runtime_error("The relative refractive index (m) should be 2-D NumPy array.");
+  if (coords.ndim() != 2)
+    throw std::runtime_error("The coordinates should be 2-D NumPy array.");
+
+
+  // allocate std::vector (to pass to the C++ function)
+  std::vector<std::vector<double> > x_cpp(x.size());
+  std::vector<std::vector<std::complex<double> > > m_cpp(m.size());
+  std::vector<double> Xc(coords.size()/3);
+  std::vector<double> Yc(coords.size()/3);
+  std::vector<double> Zc(coords.size()/3);
+
+  // copy py::array -> std::vector
+  std::memcpy(x_cpp.data(), x.data(), x.size()*sizeof(double));
+  std::memcpy(m_cpp.data(), m.data(), m.size()*sizeof(std::complex<double>));
+  std::memcpy(Xc.data(), coords.data(), coords.size()*sizeof(double)/3);
+  std::memcpy(Yc.data(), coords.data() + coords.size()*sizeof(double)/3, coords.size()*sizeof(double)/3);
+  std::memcpy(Zc.data(), coords.data() + 2*coords.size()*sizeof(double)/3, coords.size()*sizeof(double)/3);
+
+
+  // create std::vector (to get return from the C++ function)
+  std::vector<int> terms(x.shape(0));
+  std::vector<std::vector<std::vector<std::complex<double> > > > E, H;
+  E.resize(x.shape(0));
+  H.resize(x.shape(0));
+
+  for (ssize_t i = 0; i < x.shape(0); i++) {
+    E[i].resize(coords.shape(0));
+    H[i].resize(coords.shape(0));
+
+    for (int j = 0; j < 3; j++) {
+      E[i][j].resize(3);
+      H[i][j].resize(3);
+    }
+
+    terms[i] = nmie::nField(x.shape(1), pl, x_cpp[i], m_cpp[i], nmax, coords.shape(0), Xc, Yc, Zc, E[i], H[i]);
+  }
+
+  return py::make_tuple(terms, E, H);
+}
+

src/pb11_nmie.hpp

@@ -0,0 +1,54 @@
+//**********************************************************************************//
+//    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 <complex>
+#include <vector>
+
+
+namespace py = pybind11;
+
+
+py::tuple scattcoeffs(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                      py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                      int nmax, int pl);
+
+
+py::tuple scattnlay(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                    py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                    py::array_t<double, py::array::c_style | py::array::forcecast> theta,
+                    int nmax, int pl);
+
+
+py::tuple fieldnlay(py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                    py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> m,
+                    py::array_t<double, py::array::c_style | py::array::forcecast> coords,
+                    int nmax, int pl);
+

src/pb11_wrapper.cc

@@ -0,0 +1,22 @@
+#include <pybind11/pybind11.h>
+#include "pb11_nmie.hpp"
+
+namespace py = pybind11;
+
+// wrap as Python module
+PYBIND11_MODULE(scattnlay, m)
+{
+  m.doc() = "The Python version of scattnlay";
+
+  m.def("scattcoeffs", &scattcoeffs,
+        "Calculate the scattering coefficients, required to calculate both the near- and far-field parameters.",
+        py::arg("x"), py::arg("m"), py::arg("nmax"), py::arg("pl"));
+
+  m.def("scattnlay", &scattnlay,
+        "Calculate the scattering scattering parameters and amplitudes.",
+        py::arg("x"), py::arg("m"), py::arg("theta"), py::arg("nmax"), py::arg("pl"));
+
+  m.def("fieldnlay", &fieldnlay,
+        "Calculate the complex electric and magnetic field in the surroundings and inside the particle.",
+        py::arg("x"), py::arg("m"), py::arg("coords"), py::arg("nmax"), py::arg("pl"));
+}

src/pb11_wrapper_mp.cc

@@ -0,0 +1,22 @@
+#include <pybind11/pybind11.h>
+#include "pb11_nmie.hpp"
+
+namespace py = pybind11;
+
+// wrap as Python module
+PYBIND11_MODULE(scattnlay_mp, m)
+{
+  m.doc() = "The Python version of scattnlay";
+
+  m.def("scattcoeffs", &scattcoeffs,
+        "Calculate the scattering coefficients, required to calculate both the near- and far-field parameters.",
+        py::arg("x"), py::arg("m"), py::arg("nmax"), py::arg("pl"));
+
+  m.def("scattnlay", &scattnlay,
+        "Calculate the scattering scattering parameters and amplitudes.",
+        py::arg("x"), py::arg("m"), py::arg("theta"), py::arg("nmax"), py::arg("pl"));
+
+  m.def("fieldnlay", &fieldnlay,
+        "Calculate the complex electric and magnetic field in the surroundings and inside the particle.",
+        py::arg("x"), py::arg("m"), py::arg("coords"), py::arg("nmax"), py::arg("pl"));
+}