Bläddra i källkod

initial support of templated getters

Konstantin Ladutenko 3 år sedan
förälder
incheckning
6abacb078b

+ 1 - 0
examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -67,6 +67,7 @@ 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)
 
 scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
 zero = np.zeros(npts*npts, dtype = np.float64)

+ 22 - 4
scattnlay/main.py

@@ -167,17 +167,35 @@ 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_
-    return scattnlay2_(x, m, theta, nmax=nmax, pl=pl)
-        # from scattnlay_dp import mie_dp as mie_
-    # mie = mie_()
-    # a = mie.GetPECLayer()
+    from scattnlay_dp import mie_dp as mie_
+    mie = mie_()
+    mie.SetLayersSize(x)
+    mie.SetLayersIndex(m)
+    mie.SetAngles(theta)
+    mie.SetPECLayer(pl)
+    mie.SetMaxTerms(nmax)
+    mie.RunMieCalculation()
+    Qext =  mie.GetQext()
+    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)
 
 def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     """

+ 28 - 7
src/nmie-basic.hpp

@@ -64,14 +64,14 @@ namespace nmie {
   // Returns previously calculated Qext                                     //
   // ********************************************************************** //
   template <typename FloatType>
-  FloatType MultiLayerMie<FloatType>::GetQext() {
+  template <typename outputType>
+  outputType MultiLayerMie<FloatType>::GetQext() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    return Qext_;
+    return static_cast<outputType>(Qext_);
   }
 
-
-  // ********************************************************************** //
+// ********************************************************************** //
   // Returns previously calculated Qabs                                     //
   // ********************************************************************** //
   template <typename FloatType>
@@ -189,7 +189,9 @@ namespace nmie {
   }
 
 
-  // ********************************************************************** //
+
+
+// ********************************************************************** //
   // Modify refractive index of all layers                                  //
   // ********************************************************************** //
   template <typename FloatType>
@@ -198,7 +200,6 @@ namespace nmie {
     refractive_index_ = index;
   }
 
-
   // ********************************************************************** //
   // Modify coordinates for field calculation                               //
   // ********************************************************************** //
@@ -233,8 +234,28 @@ namespace nmie {
     nmax_preset_ = nmax;
   }
 
+  // Python interface
+  template <typename FloatType>
+  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);
+    SetLayersSize(ConvertVector<FloatType>(layer_size_dp));
+  }
 
-  // ********************************************************************** //
+  template <typename FloatType>
+  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);
+    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);
+   SetAngles(ConvertVector<FloatType>(angles_dp));
+ }
+
+// ********************************************************************** //
   // Get total size parameter of particle                                   //
   // ********************************************************************** //
   template <typename FloatType>

+ 0 - 6
src/nmie-pybind11.cc

@@ -84,12 +84,6 @@ py::array Vector2DComplex2Py(const std::vector<std::vector<T > > &x) {
 }
 
 
-template <typename T>
-std::vector<T> Py2Vector(const py::array_t<T> &py_x) {
-  std::vector<T> c_x(py_x.size());
-  std::memcpy(c_x.data(), py_x.data(), py_x.size()*sizeof(T));
-  return c_x;
-}
 
 
 py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,

+ 20 - 20
src/nmie.cc

@@ -190,33 +190,33 @@ namespace nmie {
     if (Theta.size() != nTheta)
         throw std::invalid_argument("Declared number of sample for Theta is not correct!");
     try {
-      MultiLayerMie<FloatType> ml_mie;
-      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
-      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
-      ml_mie.SetAngles(ConvertVector<FloatType>(Theta));
-      ml_mie.SetPECLayer(pl);
-      ml_mie.SetMaxTerms(nmax);
-      ml_mie.SetModeNmaxAndType(mode_n, mode_type);
+      MultiLayerMie<FloatType> mie;
+      mie.SetLayersSize(ConvertVector<FloatType>(x));
+      mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
+      mie.SetAngles(ConvertVector<FloatType>(Theta));
+      mie.SetPECLayer(pl);
+      mie.SetMaxTerms(nmax);
+      mie.SetModeNmaxAndType(mode_n, mode_type);
 
-      ml_mie.RunMieCalculation();
+      mie.RunMieCalculation();
 
       // std::cout
       // 	<< std::setprecision(std::numeric_limits<FloatType>::digits10)
       // 	<< "Qext = "
-      // 	<< ml_mie.GetQext()
+      // 	<< mie.GetQext()
       // 	<< std::endl;
 
-      *Qext = static_cast<double>(ml_mie.GetQext());
-      *Qsca = static_cast<double>(ml_mie.GetQsca());
-      *Qabs = static_cast<double>(ml_mie.GetQabs());
-      *Qbk = static_cast<double>(ml_mie.GetQbk());
-      *Qpr = static_cast<double>(ml_mie.GetQpr());
-      *g = static_cast<double>(ml_mie.GetAsymmetryFactor());
-      *Albedo = static_cast<double>(ml_mie.GetAlbedo());
-      S1 = ConvertComplexVector<double>(ml_mie.GetS1());
-      S2 = ConvertComplexVector<double>(ml_mie.GetS2());
+      *Qext = static_cast<double>(mie.GetQext());
+      *Qsca = static_cast<double>(mie.GetQsca());
+      *Qabs = static_cast<double>(mie.GetQabs());
+      *Qbk = static_cast<double>(mie.GetQbk());
+      *Qpr = static_cast<double>(mie.GetQpr());
+      *g = static_cast<double>(mie.GetAsymmetryFactor());
+      *Albedo = static_cast<double>(mie.GetAlbedo());
+      S1 = ConvertComplexVector<double>(mie.GetS1());
+      S2 = ConvertComplexVector<double>(mie.GetS2());
 
-      return ml_mie.GetMaxTerms();
+      return mie.GetMaxTerms();
     } catch(const std::invalid_argument &ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
@@ -415,7 +415,7 @@ namespace nmie {
       ml_mie.SetFieldCoords({ConvertVector<FloatType>(Xp_vec),
 	    ConvertVector<FloatType>(Yp_vec),
 	    ConvertVector<FloatType>(Zp_vec) });
-      ml_mie.SetMaxTerms(nmax);
+      if (nmax != -1) ml_mie.SetMaxTerms(nmax);
       ml_mie.SetModeNmaxAndType(mode_n, mode_type);
 
       ml_mie.RunFieldCalculation();

+ 17 - 3
src/nmie.hpp

@@ -37,10 +37,21 @@
 #include <cstdlib>
 #include <iostream>
 #include <vector>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
 #include "nmie-precision.hpp"
 #ifdef MULTI_PRECISION
 #include <boost/math/constants/constants.hpp>
 #endif
+namespace py = pybind11;
+
+template <typename T>
+std::vector<T> Py2Vector(const py::array_t<T> &py_x) {
+  std::vector<T> c_x(py_x.size());
+  std::memcpy(c_x.data(), py_x.data(), py_x.size()*sizeof(T));
+  return c_x;
+}
+
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl,
                   const std::vector<double> &x, const std::vector<std::complex<double> > &m,
@@ -151,7 +162,7 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
     void calcExpanCoeffs();
 
     // Return calculation results
-    FloatType GetQext();
+    template <typename outputType = FloatType> outputType GetQext();
     FloatType GetQsca();
     FloatType GetQabs();
     FloatType GetQbk();
@@ -172,12 +183,15 @@ 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);
     // 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);
     void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml, unsigned int &l);
     void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml);
-    // Modify scattering (theta) angles
-    void SetAngles(const std::vector<FloatType> &angles);
+    // 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);
     // Modify coordinates for field calculation
     void SetFieldCoords(const std::vector< std::vector<FloatType> > &coords);
     // Modify index of PEC layer

+ 27 - 4
src/pb11_wrapper.cc

@@ -1,9 +1,13 @@
 #include <string>
-#include <pybind11/pybind11.h>
 #include "nmie-pybind11.hpp"
 #include "nmie.hpp"
+#include "nmie-basic.hpp"
 
-namespace py = pybind11;
+//py::class_<Pet>(m, "Pet")
+//.def(py::init<const std::string &, int>())
+//.def("set", static_cast<void (Pet::*)(int)>(&Pet::set), "Set the pet's age")
+//.def("set", static_cast<void (Pet::*)(const std::string &)>(&Pet::set), "Set the pet's name");
+//
 
 template<typename T>
 void declare_nmie(py::module &m, std::string &typestr) {
@@ -12,8 +16,27 @@ void declare_nmie(py::module &m, std::string &typestr) {
   py::class_<mie_typed>(m, pyclass_name.c_str(), py::buffer_protocol(), py::dynamic_attr())
       .def(py::init<>())
       .def("GetPECLayer", &mie_typed::GetPECLayer)
-//      .def("width",     &mie_typed::width)
-//      .def("height",    &mie_typed::height)
+      .def("SetPECLayer", &mie_typed::SetPECLayer)
+      .def("SetMaxTerms", &mie_typed::SetMaxTerms)
+      .def("GetMaxTerms", &mie_typed::GetMaxTerms)
+      .def("SetModeNmaxAndType", &mie_typed::SetModeNmaxAndType)
+      .def("RunMieCalculation", &mie_typed::RunMieCalculation)
+      .def("SetLayersSize", static_cast<
+               void (mie_typed::*)
+                   (const py::array_t<double, py::array::c_style | py::array::forcecast>&)>
+           (&mie_typed::SetLayersSize)
+      )
+      .def("SetLayersIndex", static_cast<
+               void (mie_typed::*)
+                   (const py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> &)>
+           (&mie_typed::SetLayersIndex)
+      )
+      .def("SetAngles", static_cast<
+               void (mie_typed::*)
+                   (const py::array_t<double, py::array::c_style | py::array::forcecast>&)>
+           (&mie_typed::SetAngles)
+      )
+      .def("GetQext", &mie_typed::GetQext<double>)
       ;
 }