scattcoeffs is now emulated on Python via C++ object calls

examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -31,8 +31,7 @@
 #    along with this program.  If not, see <http:#www.gnu.org/licenses/>.
 
 import scattnlay
-from scattnlay import fieldnlay
-from scattnlay import scattnlay
+from scattnlay import fieldnlay, scattnlay, scattcoeffs, expancoeffs
 import numpy as np
 from matplotlib import pyplot as plt
 
@@ -40,7 +39,7 @@ import inspect
 print("Using scattnlay from ", inspect.getfile(scattnlay))
 
 npts = 351
-# npts = 51
+# npts = 11
 factor = 3 # plot extent compared to sphere radius
 index_H2O = 1.33+(1e-6)*1j
 
@@ -67,8 +66,12 @@ print("   Qsca = " + str(Qsca)+" terms = "+str(terms))
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
     np.array([x]), np.array([m]), mp=True)
 print("mp Qsca = " + str(Qsca)+" terms = "+str(terms))
-exit(1)
 
+terms, a,b = scattcoeffs(np.array([x]), np.array([m]))
+print(a)
+print(b)
+
+exit(1)
 scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
 zero = np.zeros(npts*npts, dtype = np.float64)
 

scattnlay/main.py

@@ -32,6 +32,25 @@
 
 import numpy as np
 
+
+def scattcoeffs_(x, m, nmax=-1, pl=-1, mp=False):
+    if mp:
+        from scattnlay_mp import mie_mp as mie_
+    else:
+        # from scattnlay_dp import mie_dp as mie_
+        from scattnlay_mp import mie_mp as mie_
+    mie = mie_()
+    mie.SetLayersSize(x)
+    mie.SetLayersIndex(m)
+    mie.SetPECLayer(pl)
+    mie.SetMaxTerms(nmax)
+    mie.calcScattCoeffs()
+    terms = mie.GetMaxTerms()
+    a = mie.GetAn()
+    b = mie.GetBn()
+    return terms, a, b
+
+
 def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     """
     scattcoeffs(x, m[, nmax, pl, mp])
@@ -53,16 +72,11 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
         an, bn: Complex scattering coefficients
     """
 
-    if mp:
-        from scattnlay_mp import scattcoeffs as scattcoeffs_
-    else:
-        from scattnlay_dp import scattcoeffs as scattcoeffs_
-
     if len(m.shape) != 1 and len(m.shape) != 2:
         raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
     if len(x.shape) == 1:
         if len(m.shape) == 1:
-            return scattcoeffs_(x, m, nmax=nmax, pl=pl)
+            return scattcoeffs_(x, m, nmax=nmax, pl=pl, mp=mp)
         else:
             raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
     elif len(x.shape) != 2:
@@ -82,7 +96,7 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     bn = np.zeros((0, nstore), dtype=complex)
 
     for i, xi in enumerate(x):
-        terms[i], a, b = scattcoeffs_(xi, m[i], nmax=nmax, pl=pl)
+        terms[i], a, b = scattcoeffs_(xi, m[i], nmax=nmax, pl=pl, mp=mp)
 
         if terms[i] > nstore:
             nstore = terms[i]
@@ -220,7 +234,7 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
         raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
     if len(x.shape) == 1:
         if len(m.shape) == 1:
-            return scattnlay_(x, m, theta, nmax=nmax, pl=pl)
+            return scattnlay_(x, m, theta, nmax=nmax, pl=pl, mp=mp)
         else:
             raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
     elif len(x.shape) != 2:
@@ -250,6 +264,25 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
 #scattnlay()
 
 
+def fieldnlay_(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
+    if mp:
+        from scattnlay_mp import mie_mp as mie_
+    else:
+        from scattnlay_dp import mie_dp as mie_
+        # from scattnlay_mp import mie_mp as mie_
+    mie = mie_()
+    mie.SetLayersSize(x)
+    mie.SetLayersIndex(m)
+    mie.SetPECLayer(pl)
+    mie.SetMaxTerms(nmax)
+    mie.SetFieldCoords(xp, yp, zp)
+    mie.RunFieldCalculation()
+    terms = mie.GetMaxTerms()
+    E = mie.GetFieldE()
+    H = mie.GetFieldH()
+    return terms, E, H
+
+
 def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     """
     fieldnlay(x, m, xp, yp, zp[, nmax, pl, mp])
@@ -278,16 +311,12 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
            (or structure) and correct it for the following ones
     """
 
-    if mp:
-        from scattnlay_mp import fieldnlay as fieldnlay_
-    else:
-        from scattnlay_dp import fieldnlay as fieldnlay_
 
     if len(m.shape) != 1 and len(m.shape) != 2:
         raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
     if len(x.shape) == 1:
         if len(m.shape) == 1:
-            return fieldnlay_(x, m, xp, yp, zp, nmax=nmax, pl=pl)
+            return fieldnlay_(x, m, xp, yp, zp, nmax=nmax, pl=pl, mp=mp)
         else:
             raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
     elif len(x.shape) != 2:
@@ -304,7 +333,7 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     for i, xi in enumerate(x):
         # (2020/05/12) We assume that the coordinates are referred to the first wavelength
         #              (or structure) and correct it for the following ones
-        terms[i], E[i], H[i] = fieldnlay_(xi, m[i], xp*xi[-1]/x[0, -1], yp*xi[-1]/x[0, -1], zp*xi[-1]/x[0, -1], nmax=nmax, pl=pl)
+        terms[i], E[i], H[i] = fieldnlay_(xi, m[i], xp*xi[-1]/x[0, -1], yp*xi[-1]/x[0, -1], zp*xi[-1]/x[0, -1], nmax=nmax, pl=pl, mp=mp)
 
     return terms, E, H
 #fieldnlay()

src/nmie-basic.hpp

@@ -163,14 +163,7 @@ namespace nmie {
       throw std::invalid_argument("You should run calculations before result request!");
     return S2_;
   }
-template <typename FloatType> template <typename outputType>
-py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS1() {
-  return VectorComplex2Py<FloatType,outputType>(GetS1());
-}
-template <typename FloatType> template <typename outputType>
-py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS2() {
-  return VectorComplex2Py<FloatType,outputType>(GetS2());
-}
+
 
 
 // ********************************************************************** //
@@ -248,30 +241,6 @@ py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS2() {
     nmax_preset_ = nmax;
   }
 
-  // Python interface
-  template <typename FloatType>
-  template <typename inputType>
-  void MultiLayerMie<FloatType>::SetLayersSize(
-      const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_layer_size) {
-    auto layer_size_dp = Py2Vector<inputType>(py_layer_size);
-    SetLayersSize(ConvertVector<FloatType>(layer_size_dp));
-  }
-
-  template <typename FloatType>
-  template <typename inputType>
-  void MultiLayerMie<FloatType>::SetLayersIndex(
-      const py::array_t<std::complex<inputType>, py::array::c_style | py::array::forcecast> &py_index) {
-    auto index_dp = Py2Vector<std::complex<inputType> >(py_index);
-    SetLayersIndex(ConvertComplexVector<FloatType>(index_dp));
-  }
-
- template <typename FloatType>
- template <typename inputType>
- void MultiLayerMie<FloatType>::SetAngles(const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_angles) {
-   auto angles_dp = Py2Vector<inputType>(py_angles);
-   SetAngles(ConvertVector<FloatType>(angles_dp));
- }
-
 // ********************************************************************** //
   // Get total size parameter of particle                                   //
   // ********************************************************************** //
@@ -489,10 +458,10 @@ py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS2() {
   void MultiLayerMie<FloatType>::calcPiTauAllTheta(const double from_Theta, const double to_Theta,
                                                    std::vector<std::vector<FloatType> > &Pi,
                                                    std::vector<std::vector<FloatType> > &Tau) {
-    auto perimeter_points = Pi.size();
+    const unsigned int perimeter_points = Pi.size();
     for (auto &val:Pi) val.resize(available_maximal_nmax_);
     for (auto &val:Tau) val.resize(available_maximal_nmax_);
-    double delta_Theta = eval_delta<FloatType>(perimeter_points, from_Theta, to_Theta);
+    double delta_Theta = eval_delta<double>(perimeter_points, from_Theta, to_Theta);
     for (unsigned int i=0; i < perimeter_points; i++) {
       auto Theta = static_cast<FloatType>(from_Theta + i*delta_Theta);
       // Calculate angular functions Pi and Tau
@@ -901,5 +870,73 @@ py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS2() {
     isMieCalculated_ = true;
   }
 
+  // Python interface
+  template <typename FloatType> template <typename inputType>
+  void MultiLayerMie<FloatType>::SetLayersSize(
+      const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_layer_size) {
+    auto layer_size_dp = Py2Vector<inputType>(py_layer_size);
+    SetLayersSize(ConvertVector<FloatType>(layer_size_dp));
+  }
+
+  template <typename FloatType> template <typename inputType>
+  void MultiLayerMie<FloatType>::SetLayersIndex(
+      const py::array_t<std::complex<inputType>, py::array::c_style | py::array::forcecast> &py_index) {
+    auto index_dp = Py2Vector<std::complex<inputType> >(py_index);
+    SetLayersIndex(ConvertComplexVector<FloatType>(index_dp));
+  }
+
+  template <typename FloatType> template <typename inputType>
+  void MultiLayerMie<FloatType>::SetAngles(const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_angles) {
+    auto angles_dp = Py2Vector<inputType>(py_angles);
+    SetAngles(ConvertVector<FloatType>(angles_dp));
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS1() {
+    return VectorComplex2Py<FloatType,outputType>(GetS1());
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetS2() {
+    return VectorComplex2Py<FloatType,outputType>(GetS2());
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetAn() {
+    return VectorComplex2Py<FloatType,outputType>(GetAn());
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array_t< std::complex<outputType>>  MultiLayerMie<FloatType>::GetBn() {
+    return VectorComplex2Py<FloatType,outputType>(GetBn());
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array MultiLayerMie<FloatType>::GetFieldE() {
+    return Vector2DComplex2Py<std::complex<outputType> >(
+        ConvertComplexVectorVector<outputType>(GetFieldE())
+    );
+  }
+
+  template <typename FloatType> template <typename outputType>
+  py::array MultiLayerMie<FloatType>::GetFieldH() {
+    return Vector2DComplex2Py<std::complex<outputType> >(
+        ConvertComplexVectorVector<outputType>(GetFieldH())
+    );
+  }
+
+  template <typename FloatType>
+  void MultiLayerMie<FloatType>::SetFieldCoords(
+      const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Xp,
+      const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Yp,
+      const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Zp) {
+    auto c_Xp = Py2Vector<double>(py_Xp);
+    auto c_Yp = Py2Vector<double>(py_Yp);
+    auto c_Zp = Py2Vector<double>(py_Zp);
+    SetFieldCoords({ConvertVector<FloatType>(c_Xp),
+                    ConvertVector<FloatType>(c_Yp),
+                    ConvertVector<FloatType>(c_Zp) });
+  }
+
 }  // end of namespace nmie
 #endif  // SRC_NMIE_BASIC_HPP_

src/nmie-nearfield.hpp

@@ -408,9 +408,8 @@ namespace nmie {
     convertFieldsFromSphericalToCartesian();
   }  //  end of MultiLayerMie::RunFieldCalculation()
 
-
 template <typename FloatType>
-double eval_delta(const int steps, const double from_value, const double to_value) {
+double eval_delta(const unsigned int steps, const double from_value, const double to_value) {
   auto delta = std::abs(from_value - to_value);
   if (steps < 2) return delta;
   delta /= static_cast<double>(steps-1);
@@ -461,15 +460,15 @@ void MultiLayerMie<FloatType>::calcMieSeriesNeededToConverge(const FloatType Rho
 
 
 template <typename FloatType>
-void MultiLayerMie<FloatType>::calcRadialOnlyDependantFunctions(const FloatType from_Rho, const FloatType to_Rho,
+void MultiLayerMie<FloatType>::calcRadialOnlyDependantFunctions(const double from_Rho, const double to_Rho,
                                                                 const bool isIgnoreAvailableNmax,
                                                                 std::vector<std::vector<std::complex<FloatType> > > &Psi,
                                                                 std::vector<std::vector<std::complex<FloatType> > > &D1n,
                                                                 std::vector<std::vector<std::complex<FloatType> > > &Zeta,
                                                                 std::vector<std::vector<std::complex<FloatType> > > &D3n) {
-  unsigned int radius_points = Psi.size();
+  auto radius_points = Psi.size();
   std::vector<std::vector<std::complex<FloatType> > > PsiZeta(radius_points);
-  double delta_Rho = eval_delta<FloatType>(radius_points, from_Rho, to_Rho);
+  double delta_Rho = eval_delta<double>(radius_points, from_Rho, to_Rho);
   for (unsigned int j=0; j < radius_points; j++) {
     auto Rho = static_cast<FloatType>(from_Rho + j*delta_Rho);
 //    if (Rho < 1e-5) Rho = 1e-5; // TODO do we need this?.
@@ -518,9 +517,9 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int outer_arc_poin
   calcRadialOnlyDependantFunctions(from_Rho, to_Rho, isIgnoreAvailableNmax,
                                    Psi, D1n, Zeta, D3n);
 
-  double delta_Rho = eval_delta<FloatType>(radius_points, from_Rho, to_Rho);
-  double delta_Theta = eval_delta<FloatType>(outer_arc_points, from_Theta, to_Theta);
-  double delta_Phi = eval_delta<FloatType>(radius_points, from_Phi, to_Phi);
+  double delta_Rho = eval_delta<double>(radius_points, from_Rho, to_Rho);
+  double delta_Theta = eval_delta<double>(outer_arc_points, from_Theta, to_Theta);
+  double delta_Phi = eval_delta<double>(radius_points, from_Phi, to_Phi);
   Es_.clear(); Hs_.clear(); coords_polar_.clear();
   for (int j=0; j < radius_points; j++) {
     auto Rho = static_cast<FloatType>(from_Rho + j * delta_Rho);
@@ -542,5 +541,9 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int outer_arc_poin
   convertFieldsFromSphericalToCartesian();
 
 }
+
+// Python interface
+
+
 }  // end of namespace nmie
 #endif  // SRC_NMIE_NEARFIELD_HPP_

src/nmie-pybind11.cc

@@ -40,44 +40,6 @@
 namespace py = pybind11;
 
 
-// https://stackoverflow.com/questions/17294629/merging-flattening-sub-vectors-into-a-single-vector-c-converting-2d-to-1d
-template <typename T>
-std::vector<T> flatten(const std::vector<std::vector<T>> &v) {
-    std::size_t total_size = 0;
-    for (const auto &sub : v)
-        total_size += sub.size(); // I wish there was a transform_accumulate
-    std::vector<T> result;
-    result.reserve(total_size);
-    for (const auto &sub : v)
-        result.insert(result.end(), sub.begin(), sub.end());
-    return result;
-}
-
-
-template <typename T>
-py::array Vector2DComplex2Py(const std::vector<std::vector<T > > &x) {
-  size_t dim1 = x.size();
-  size_t dim2 = x[0].size();
-  auto result = flatten(x);
-  // https://github.com/tdegeus/pybind11_examples/blob/master/04_numpy-2D_cpp-vector/example.cpp
-  size_t              ndim    = 2;
-  std::vector<size_t> shape   = {dim1, dim2};
-  std::vector<size_t> strides = {sizeof(T)*dim2, sizeof(T)};
-
-  // return 2-D NumPy array
-  return py::array(py::buffer_info(
-    result.data(),                       /* data as contiguous array  */
-    sizeof(T),                           /* size of one scalar        */
-    py::format_descriptor<T>::format(),  /* data type                 */
-    ndim,                                /* number of dimensions      */
-    shape,                               /* shape of the matrix       */
-    strides                              /* strides for each axis     */
-  ));
-}
-
-
-
-
 py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
                          const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
                          const int nmax=-1, const int pl=-1) {
@@ -110,30 +72,3 @@ py::tuple py_ExpanCoeffs(const py::array_t<double, py::array::c_style | py::arra
 
   return py::make_tuple(terms, Vector2DComplex2Py(c_an), Vector2DComplex2Py(c_bn), Vector2DComplex2Py(c_cn), Vector2DComplex2Py(c_dn));
 }
-
-
-py::tuple py_fieldnlay(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,
-                       const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_m,
-                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Xp,
-                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Yp,
-                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Zp,
-                       const int nmax=-1, const int pl=-1) {
-
-  auto c_x = Py2Vector<double>(py_x);
-  auto c_m = Py2Vector< std::complex<double> >(py_m);
-  auto c_Xp = Py2Vector<double>(py_Xp);
-  auto c_Yp = Py2Vector<double>(py_Yp);
-  auto c_Zp = Py2Vector<double>(py_Zp);
-  unsigned int ncoord = py_Xp.size();
-  std::vector<std::vector<std::complex<double> > > E(ncoord);
-  std::vector<std::vector<std::complex<double> > > H(ncoord);
-  for (auto &f : E) f.resize(3);
-  for (auto &f : H) f.resize(3);
-  int L = py_x.size(), terms;
-  terms = nmie::nField(L, pl, c_x, c_m, nmax, nmie::Modes::kAll, nmie::Modes::kAll, ncoord, c_Xp, c_Yp, c_Zp, E, H);
-  auto py_E = Vector2DComplex2Py<std::complex<double> >(E);
-  auto py_H = Vector2DComplex2Py<std::complex<double> >(H);
-  return py::make_tuple(terms, py_E, py_H);
-}
-
-

src/nmie.cc

@@ -415,7 +415,7 @@ namespace nmie {
       ml_mie.SetFieldCoords({ConvertVector<FloatType>(Xp_vec),
 	    ConvertVector<FloatType>(Yp_vec),
 	    ConvertVector<FloatType>(Zp_vec) });
-      if (nmax != -1) ml_mie.SetMaxTerms(nmax);
+      ml_mie.SetMaxTerms(nmax);
       ml_mie.SetModeNmaxAndType(mode_n, mode_type);
 
       ml_mie.RunFieldCalculation();

src/nmie.hpp

@@ -37,8 +37,11 @@
 #include <cstdlib>
 #include <iostream>
 #include <vector>
+
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
+#include <pybind11/complex.h>
+
 #include "nmie-precision.hpp"
 #ifdef MULTI_PRECISION
 #include <boost/math/constants/constants.hpp>
@@ -53,8 +56,8 @@ std::vector<T> Py2Vector(const py::array_t<T> &py_x) {
 }
 
 template <typename inputType=double, typename outputType=double>
-py::array_t< std::complex<outputType>> VectorComplex2Py(const std::vector<std::complex<inputType> > &c_f) {
-  auto c_x = nmie::ConvertComplexVector<outputType, inputType>(c_f);
+py::array_t< std::complex<outputType>> VectorComplex2Py(const std::vector<std::complex<inputType> > &cf_x) {
+  auto c_x = nmie::ConvertComplexVector<outputType, inputType>(cf_x);
   auto py_x = py::array_t< std::complex<outputType>>(c_x.size());
   auto py_x_buffer = py_x.request();
   auto *py_x_ptr = (std::complex<outputType> *) py_x_buffer.ptr;
@@ -62,6 +65,41 @@ py::array_t< std::complex<outputType>> VectorComplex2Py(const std::vector<std::c
   return py_x;
 }
 
+// https://stackoverflow.com/questions/17294629/merging-flattening-sub-vectors-into-a-single-vector-c-converting-2d-to-1d
+template <typename T>
+std::vector<T> flatten(const std::vector<std::vector<T>> &v) {
+  std::size_t total_size = 0;
+  for (const auto &sub : v)
+    total_size += sub.size(); // I wish there was a transform_accumulate
+  std::vector<T> result;
+  result.reserve(total_size);
+  for (const auto &sub : v)
+    result.insert(result.end(), sub.begin(), sub.end());
+  return result;
+}
+
+
+template <typename T>
+py::array Vector2DComplex2Py(const std::vector<std::vector<T > > &x) {
+  size_t dim1 = x.size();
+  size_t dim2 = x[0].size();
+  auto result = flatten(x);
+  // https://github.com/tdegeus/pybind11_examples/blob/master/04_numpy-2D_cpp-vector/example.cpp
+  size_t              ndim    = 2;
+  std::vector<size_t> shape   = {dim1, dim2};
+  std::vector<size_t> strides = {sizeof(T)*dim2, sizeof(T)};
+
+  // return 2-D NumPy array
+  return py::array(py::buffer_info(
+      result.data(),                       /* data as contiguous array  */
+      sizeof(T),                           /* size of one scalar        */
+      py::format_descriptor<T>::format(),  /* data type                 */
+      ndim,                                /* number of dimensions      */
+      shape,                               /* shape of the matrix       */
+      strides                              /* strides for each axis     */
+  ));
+}
+
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl,
                   const std::vector<double> &x, const std::vector<std::complex<double> > &m,
@@ -79,7 +117,7 @@ namespace nmie {
 
 //helper functions
 template <typename FloatType>
-double eval_delta(const int steps, const double from_value, const double to_value);
+double eval_delta(const unsigned int steps, const double from_value, const double to_value);
 
 
 template<class T> inline T pow2(const T value) {return value*value;}
@@ -186,6 +224,8 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
 
     std::vector<std::complex<FloatType> > GetAn(){return an_;};
     std::vector<std::complex<FloatType> > GetBn(){return bn_;};
+    template <typename outputType> py::array_t< std::complex<outputType>>  GetAn();
+    template <typename outputType> py::array_t< std::complex<outputType>>  GetBn();
 
     std::vector< std::vector<std::complex<FloatType> > > GetLayerAn(){return aln_;};
     std::vector< std::vector<std::complex<FloatType> > > GetLayerBn(){return bln_;};
@@ -209,6 +249,9 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
     void SetAngles(const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_angles);
     // Modify coordinates for field calculation
     void SetFieldCoords(const std::vector< std::vector<FloatType> > &coords);
+    void SetFieldCoords(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Xp,
+                        const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Yp,
+                        const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Zp);
     // Modify index of PEC layer
     void SetPECLayer(int layer_position = 0);
     // Modify the mode taking into account for evaluation of output variables
@@ -240,6 +283,9 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
 
     std::vector<std::vector< std::complex<FloatType> > > GetFieldE(){return E_;};   // {X[], Y[], Z[]}
     std::vector<std::vector< std::complex<FloatType> > > GetFieldH(){return H_;};
+    template <typename outputType> py::array GetFieldE();
+    template <typename outputType> py::array GetFieldH();
+
     // Get fields in spherical coordinates.
     std::vector<std::vector< std::complex<FloatType> > > GetFieldEs(){return E_;};   // {rho[], teha[], phi[]}
     std::vector<std::vector< std::complex<FloatType> > > GetFieldHs(){return H_;};
@@ -320,8 +366,8 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
                            const double to_Theta,
                            std::vector<std::vector<FloatType>> &Pi,
                            std::vector<std::vector<FloatType>> &Tau);
-    void calcRadialOnlyDependantFunctions(const FloatType from_Rho,
-                                          const FloatType to_Rho,
+    void calcRadialOnlyDependantFunctions(const double from_Rho,
+                                          const double to_Rho,
                                           const bool isIgnoreAvailableNmax,
                                           std::vector<std::vector<std::complex<FloatType>>> &Psi,
                                           std::vector<std::vector<std::complex<FloatType>>> &D1n,

src/pb11_wrapper.cc

@@ -2,6 +2,7 @@
 #include "nmie-pybind11.hpp"
 #include "nmie.hpp"
 #include "nmie-basic.hpp"
+#include "nmie-nearfield.hpp"
 
 //py::class_<Pet>(m, "Pet")
 //.def(py::init<const std::string &, int>())
@@ -21,6 +22,12 @@ void declare_nmie(py::module &m, std::string &typestr) {
       .def("GetMaxTerms", &mie_typed::GetMaxTerms)
       .def("SetModeNmaxAndType", &mie_typed::SetModeNmaxAndType)
       .def("RunMieCalculation", &mie_typed::RunMieCalculation)
+      .def("calcScattCoeffs", &mie_typed::calcScattCoeffs)
+      .def("calcExpanCoeffs", &mie_typed::calcExpanCoeffs)
+      .def("RunFieldCalculation", &mie_typed::RunFieldCalculation)
+      .def("RunFieldCalculationPolar", &mie_typed::RunFieldCalculationPolar)
+      .def("GetPECLayer", &mie_typed::GetPECLayer)
+
       .def("SetLayersSize", static_cast<
                void (mie_typed::*)
                    (const py::array_t<double, py::array::c_style | py::array::forcecast>&)>
@@ -36,6 +43,15 @@ void declare_nmie(py::module &m, std::string &typestr) {
                    (const py::array_t<double, py::array::c_style | py::array::forcecast>&)>
            (&mie_typed::SetAngles)
       )
+      .def("SetFieldCoords", static_cast<
+               void (mie_typed::*)
+                   (
+                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Xp,
+                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Yp,
+                       const py::array_t<double, py::array::c_style | py::array::forcecast> &py_Zp
+                   )>
+           (&mie_typed::SetFieldCoords)
+      )
       .def("GetQext", &mie_typed::GetQext<double>)
       .def("GetQsca", &mie_typed::GetQsca<double>)
       .def("GetQabs", &mie_typed::GetQabs<double>)
@@ -45,6 +61,10 @@ void declare_nmie(py::module &m, std::string &typestr) {
       .def("GetAlbedo", &mie_typed::GetAlbedo<double>)
       .def("GetS1", &mie_typed::GetS1<double>)
       .def("GetS2", &mie_typed::GetS2<double>)
+      .def("GetAn", &mie_typed::GetAn<double>)
+      .def("GetBn", &mie_typed::GetBn<double>)
+      .def("GetFieldE", &mie_typed::GetFieldE<double>)
+      .def("GetFieldH", &mie_typed::GetFieldH<double>)
       ;
 }
 
@@ -60,15 +80,8 @@ PYBIND11_MODULE(scattnlay_dp, m)
   m.doc() = "The Python version of scattnlay";
   declare_nmie<nmie::FloatType>(m, precision_name);
 
-  m.def("scattcoeffs", &py_ScattCoeffs,
-        "Calculate the scattering coefficients, required to calculate both the near- and far-field parameters.",
-        py::arg("x"), py::arg("m"), py::arg("nmax") = -1, py::arg("pl") = -1);
-
   m.def("expancoeffs", &py_ExpanCoeffs,
         "Calculate the expansion coefficients, required to calculate the near-field parameters.",
         py::arg("x"), py::arg("m"), py::arg("nmax") = -1, py::arg("pl") = -1);
-
-    m.def("fieldnlay", &py_fieldnlay,
-        "Calculate the complex electric and magnetic field in the surroundings and inside the particle.",
-        py::arg("x"), py::arg("m"), py::arg("xp"), py::arg("yp"), py::arg("zp"), py::arg("nmax") = -1, py::arg("pl") = -1);
 }
+