remove py_scattnlay from C++ (now it is emulated in Python via class bindings)

scattnlay/main.py

@@ -166,21 +166,9 @@ def expancoeffs(x, m, nmax=-1, pl=-1, mp=False):
 
 def scattnlay_(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     if mp:
-        from scattnlay_mp import scattnlay as scattnlay2_
         from scattnlay_mp import mie_mp as mie_
-        mie = mie_()
-        mie.SetLayersSize(x)
-        mie.SetLayersIndex(m)
-        mie.SetAngles(theta)
-        mie.SetPECLayer(pl)
-        mie.SetMaxTerms(nmax)
-        mie.RunMieCalculation()
-        # Qext = mie.GetQext()
-        # print(Qext)
-        return scattnlay2_(x, m, theta, nmax=nmax, pl=pl)
-
-    from scattnlay_dp import scattnlay as scattnlay2_
-    from scattnlay_dp import mie_dp as mie_
+    else:
+        from scattnlay_dp import mie_dp as mie_
     mie = mie_()
     mie.SetLayersSize(x)
     mie.SetLayersIndex(m)
@@ -189,13 +177,16 @@ def scattnlay_(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     mie.SetMaxTerms(nmax)
     mie.RunMieCalculation()
     Qext =  mie.GetQext()
+    Qsca =  mie.GetQsca()
+    Qabs =  mie.GetQabs()
+    Qbk =  mie.GetQbk()
+    Qpr =  mie.GetQpr()
+    g =  mie.GetAsymmetryFactor()
+    Albedo =  mie.GetAlbedo()
     terms = mie.GetMaxTerms()
-    print(Qext)
-
-    # terms[i], Qext[i], Qsca[i], Qabs[i], Qbk[i], Qpr[i], g[i], Albedo[i], S1[i], S2[i] = scattnlay_(xi, m[i], theta, nmax=nmax, pl=pl)
-    #
-    # return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
-    return scattnlay2_(x, m, theta, nmax=nmax, pl=pl)
+    S1 = mie.GetS1()
+    S2 = mie.GetS2()
+    return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
 
 def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     """

src/nmie-basic.hpp

@@ -71,14 +71,15 @@ namespace nmie {
     return static_cast<outputType>(Qext_);
   }
 
-// ********************************************************************** //
+  // ********************************************************************** //
   // Returns previously calculated Qabs                                     //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetQabs() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetQabs() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return Qabs_;
+    return static_cast<outputType>(Qabs_);
   }
 
 
@@ -86,10 +87,11 @@ namespace nmie {
   // Returns previously calculated Qsca                                     //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetQsca() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetQsca() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return Qsca_;
+    return static_cast<outputType>(Qsca_);
   }
 
 
@@ -97,10 +99,11 @@ namespace nmie {
   // Returns previously calculated Qbk                                      //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetQbk() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetQbk() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return Qbk_;
+    return static_cast<outputType>(Qbk_);
   }
 
 
@@ -108,10 +111,11 @@ namespace nmie {
   // Returns previously calculated Qpr                                      //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetQpr() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetQpr() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return Qpr_;
+    return static_cast<outputType>(Qpr_);
   }
 
 
@@ -119,10 +123,11 @@ namespace nmie {
   // Returns previously calculated assymetry factor                         //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetAsymmetryFactor() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetAsymmetryFactor() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return asymmetry_factor_;
+    return static_cast<outputType>(asymmetry_factor_);
   }
 
 
@@ -130,10 +135,11 @@ namespace nmie {
   // Returns previously calculated Albedo                                   //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetAlbedo() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetAlbedo() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return albedo_;
+    return static_cast<outputType>(albedo_);
   }
 
 
@@ -157,9 +163,17 @@ 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());
+}
 
 
-  // ********************************************************************** //
+// ********************************************************************** //
   // Modify scattering (theta) angles                                       //
   // ********************************************************************** //
   template <typename FloatType>
@@ -236,22 +250,25 @@ namespace nmie {
 
   // Python interface
   template <typename FloatType>
+  template <typename inputType>
   void MultiLayerMie<FloatType>::SetLayersSize(
-      const py::array_t<double, py::array::c_style | py::array::forcecast> &py_layer_size) {
-    auto layer_size_dp = Py2Vector<double>(py_layer_size);
+      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<double>, py::array::c_style | py::array::forcecast> &py_index) {
-    auto index_dp = Py2Vector<std::complex<double> >(py_index);
+      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>
- void MultiLayerMie<FloatType>::SetAngles(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_angles) {
-   auto angles_dp = Py2Vector<double>(py_angles);
+ 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));
  }
 

src/nmie-pybind11.cc

@@ -39,14 +39,6 @@
 
 namespace py = pybind11;
 
-py::array_t< std::complex<double>> VectorComplex2Py(const std::vector<std::complex<double> > &c_x) {
-  auto py_x = py::array_t< std::complex<double>>(c_x.size());
-  auto py_x_buffer = py_x.request();
-  auto *py_x_ptr = (std::complex<double> *) py_x_buffer.ptr;
-  std::memcpy(py_x_ptr, c_x.data(), c_x.size()*sizeof(std::complex<double>));
-  return py_x;
-}
-
 
 // https://stackoverflow.com/questions/17294629/merging-flattening-sub-vectors-into-a-single-vector-c-converting-2d-to-1d
 template <typename T>
@@ -120,27 +112,6 @@ py::tuple py_ExpanCoeffs(const py::array_t<double, py::array::c_style | py::arra
 }
 
 
-py::tuple py_scattnlay(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_theta,
-                       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_theta = Py2Vector<double>(py_theta);
-
-  int L = py_x.size(), nTheta = c_theta.size(), terms;
-  double Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo;
-  std::vector<std::complex<double> > c_S1, c_S2;
-
-
-  terms = nmie::nMie(L, pl, c_x, c_m, nTheta, c_theta, nmax, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, c_S1, c_S2);
-
-  return py::make_tuple(terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo,
-                        VectorComplex2Py(c_S1), VectorComplex2Py(c_S2));
-}
-
-
 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,

src/nmie.hpp

@@ -52,6 +52,16 @@ std::vector<T> Py2Vector(const py::array_t<T> &py_x) {
   return c_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);
+  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;
+  std::memcpy(py_x_ptr, c_x.data(), c_x.size()*sizeof(std::complex<outputType>));
+  return py_x;
+}
+
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl,
                   const std::vector<double> &x, const std::vector<std::complex<double> > &m,
@@ -163,14 +173,16 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
 
     // Return calculation results
     template <typename outputType = FloatType> outputType GetQext();
-    FloatType GetQsca();
-    FloatType GetQabs();
-    FloatType GetQbk();
-    FloatType GetQpr();
-    FloatType GetAsymmetryFactor();
-    FloatType GetAlbedo();
+    template <typename outputType = FloatType> outputType  GetQsca();
+    template <typename outputType = FloatType> outputType  GetQabs();
+    template <typename outputType = FloatType> outputType  GetQbk();
+    template <typename outputType = FloatType> outputType  GetQpr();
+    template <typename outputType = FloatType> outputType  GetAsymmetryFactor();
+    template <typename outputType = FloatType> outputType  GetAlbedo();
     std::vector<std::complex<FloatType> > GetS1();
     std::vector<std::complex<FloatType> > GetS2();
+    template <typename outputType> py::array_t< std::complex<outputType>>  GetS1();
+    template <typename outputType> py::array_t< std::complex<outputType>>  GetS2();
 
     std::vector<std::complex<FloatType> > GetAn(){return an_;};
     std::vector<std::complex<FloatType> > GetBn(){return bn_;};
@@ -183,15 +195,18 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
     // Problem definition
     // Modify size of all layers
     void SetLayersSize(const std::vector<FloatType> &layer_size);
-    void SetLayersSize(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_layer_size);
+    template <typename inputType>
+    void SetLayersSize(const py::array_t<inputType, py::array::c_style | py::array::forcecast> &py_layer_size);
     // Modify refractive index of all layers
     void SetLayersIndex(const std::vector< std::complex<FloatType> > &index);
-    void SetLayersIndex(const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &py_index);
+    template <typename inputType>
+    void SetLayersIndex(const py::array_t<std::complex<inputType>, py::array::c_style | py::array::forcecast> &py_index);
     void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml, unsigned int &l);
     void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml);
     // Modify scattering (theta) py_angles
     void SetAngles(const std::vector<FloatType> &py_angles);
-    void SetAngles(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_angles);
+    template <typename inputType>
+    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);
     // Modify index of PEC layer

src/pb11_wrapper.cc

@@ -37,6 +37,14 @@ void declare_nmie(py::module &m, std::string &typestr) {
            (&mie_typed::SetAngles)
       )
       .def("GetQext", &mie_typed::GetQext<double>)
+      .def("GetQsca", &mie_typed::GetQsca<double>)
+      .def("GetQabs", &mie_typed::GetQabs<double>)
+      .def("GetQpr", &mie_typed::GetQpr<double>)
+      .def("GetQbk", &mie_typed::GetQbk<double>)
+      .def("GetAsymmetryFactor", &mie_typed::GetAsymmetryFactor<double>)
+      .def("GetAlbedo", &mie_typed::GetAlbedo<double>)
+      .def("GetS1", &mie_typed::GetS1<double>)
+      .def("GetS2", &mie_typed::GetS2<double>)
       ;
 }
 
@@ -60,11 +68,7 @@ PYBIND11_MODULE(scattnlay_dp, m)
         "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("scattnlay", &py_scattnlay,
-        "Calculate the scattering parameters and amplitudes.",
-        py::arg("x"), py::arg("m"), py::arg("theta") = py::array_t<double>(0), py::arg("nmax") = -1, py::arg("pl") = -1);
-
-  m.def("fieldnlay", &py_fieldnlay,
+    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);
 }