mesomie is now available from python

scattnlay/__init__.py

@@ -1,8 +1,8 @@
 #!/usr/bin/env python
 # -*- coding: UTF-8 -*-
 #
-#    Copyright (C) 2009-2019 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
-#    Copyright (C) 2013-2019 Konstantin Ladutenko <kostyfisik@gmail.com>
+#    Copyright (C) 2009-2022 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2022 Konstantin Ladutenko <kostyfisik@gmail.com>
 #
 #    This file is part of scattnlay
 #
@@ -30,4 +30,6 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay, mie, mie_mp
+from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay
+from scattnlay.main import mie, mie_mp
+from scattnlay.main import mesomie  # , mesomie_mp

scattnlay/main.py

@@ -29,18 +29,22 @@
 #
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-from scattnlay_dp import mie_dp
+from scattnlay_dp import mie_dp, mesomie_dp
 import numpy as np
 import sys
 
 mie_mp = None
+# mesomie_mp = None
 try:
     from scattnlay_mp import mie_mp as mie_mp_
     mie_mp = mie_mp_()
+    # from scattnlay_mp import mesomie_mp as mesomie_mp_
+    # mesomie_mp = mesomie_mp_()
 except:
     pass
 
 mie = mie_dp()
+mesomie = mesomie_dp()
 
 
 def scattcoeffs_(x, m, nmax=-1, pl=-1, mp=False):

src/nmie.hpp

@@ -296,8 +296,12 @@ class MultiLayerMie {
   std::vector<std::complex<FloatType>> GetS1();
   std::vector<std::complex<FloatType>> GetS2();
 
-  std::vector<std::complex<FloatType>> GetAn() { return an_; };
-  std::vector<std::complex<FloatType>> GetBn() { return bn_; };
+  std::vector<std::complex<FloatType>> GetAn() {
+    return an_;
+  };
+  std::vector<std::complex<FloatType>> GetBn() {
+    return bn_;
+  };
 
   std::vector<std::vector<std::complex<FloatType>>> GetLayerAn() {
     return aln_;
@@ -347,9 +351,13 @@ class MultiLayerMie {
   void SetMaxTerms(int nmax);
 
   // Get maximum number of terms
-  int GetMaxTerms() { return nmax_; };
+  int GetMaxTerms() {
+    return nmax_;
+  };
 
-  bool isMieCalculated() { return isMieCalculated_; };
+  bool isMieCalculated() {
+    return isMieCalculated_;
+  };
 
   // Clear layer information
   void ClearLayers();
@@ -360,25 +368,39 @@ class MultiLayerMie {
   // Get total size parameter of particle
   FloatType GetSizeParameter();
   // Returns size of all layers
-  std::vector<FloatType> GetLayersSize() { return size_param_; };
+  std::vector<FloatType> GetLayersSize() {
+    return size_param_;
+  };
   // Returns refractive index of all layers
   std::vector<std::complex<FloatType>> GetLayersIndex() {
     return refractive_index_;
   };
   // Returns scattering (theta) angles
-  std::vector<FloatType> GetAngles() { return theta_; };
+  std::vector<FloatType> GetAngles() {
+    return theta_;
+  };
   // Returns coordinates used for field calculation
-  std::vector<std::vector<FloatType>> GetFieldCoords() { return coords_; };
+  std::vector<std::vector<FloatType>> GetFieldCoords() {
+    return coords_;
+  };
   // Returns index of PEC layer
-  int GetPECLayer() { return PEC_layer_position_; };
+  int GetPECLayer() {
+    return PEC_layer_position_;
+  };
 
   std::vector<std::vector<std::complex<FloatType>>> GetFieldE() {
     return E_;
   };  // {X[], Y[], Z[]}
-  std::vector<std::vector<std::complex<FloatType>>> GetFieldH() { return H_; };
+  std::vector<std::vector<std::complex<FloatType>>> GetFieldH() {
+    return H_;
+  };
 
-  std::vector<FloatType> GetFieldEabs() { return Eabs_; };  // {X[], Y[], Z[]}
-  std::vector<FloatType> GetFieldHabs() { return Habs_; };
+  std::vector<FloatType> GetFieldEabs() {
+    return Eabs_;
+  };  // {X[], Y[], Z[]}
+  std::vector<FloatType> GetFieldHabs() {
+    return Habs_;
+  };
   bool GetFieldConvergence();
 
   // Get fields in spherical coordinates.
@@ -540,8 +562,17 @@ class MesoMie {
   std::vector<std::complex<FloatType>> GetBn() { return bn_; };
 
   FloatType Qsca_ = 0.0, Qext_ = 0.0;
-  FloatType GetQsca() { return Qsca_; };
-  FloatType GetQext() { return Qext_; };
+
+  template <typename outputType = FloatType>
+  outputType GetQsca() {
+    return static_cast<outputType>(Qsca_);
+  };
+
+  template <typename outputType = FloatType>
+  outputType GetQext() {
+    return static_cast<outputType>(Qext_);
+  };
+
   void calc_ab(FloatType R,
                std::complex<FloatType> xd,
                std::complex<FloatType> xm,

src/pb11-wrapper.cc

@@ -30,6 +30,7 @@
 //******************************************************************************
 
 #include <string>
+#include "mesomie.hpp"
 #include "nmie-basic.hpp"
 #include "nmie-nearfield.hpp"
 #include "nmie.hpp"
@@ -41,6 +42,7 @@
 
 namespace py = pybind11;
 
+//******************************************************************************
 template <typename T>
 void declare_nmie(py::module& m, const std::string& typestr) {
   using mie_typed = nmie::PyMultiLayerMie<nmie::FloatType>;
@@ -110,8 +112,22 @@ void declare_nmie(py::module& m, const std::string& typestr) {
       .def("GetLayerBn", &mie_typed::GetLayerBn<double>)
       .def("GetLayerCn", &mie_typed::GetLayerCn<double>)
       .def("GetLayerDn", &mie_typed::GetLayerDn<double>);
-}
 
+  using mesomie = nmie::MesoMie<nmie::FloatType>;
+  pyclass_name = std::string("mesomie") + typestr;
+  py::class_<mesomie>(m, pyclass_name.c_str(), py::buffer_protocol(),
+                      py::dynamic_attr())
+      .def(py::init<>())
+      .def("calc_Q", &mesomie::calc_Q)
+      .def("calc_ab", &mesomie::calc_ab, py::arg("R") = 1, py::arg("xd") = 1,
+           py::arg("xm") = 1, py::arg("eps_d") = 1, py::arg("eps_m") = 1,
+           py::arg("d_parallel") = 0, py::arg("d_perp") = 0)
+      .def("GetQext", &mesomie::GetQext<double>)
+      .def("GetQsca", &mesomie::GetQsca<double>);
+
+}  // end of declare_nmie(..)
+
+//******************************************************************************
 // wrap as Python module
 #ifdef MULTI_PRECISION
 std::string precision_name = "_mp";

tests/test_bulk_sphere.cc

@@ -57,8 +57,8 @@ std::vector<std::tuple<double, std::complex<double>, double, double, char> >
         {10000, {10, 10}, 2.005914, 1.795393, 'm'},
     };
 //******************************************************************************
-TEST(BulkSphere, DISABLED_HandlesInput) {
-  // TEST(BulkSphere, HandlesInput) {
+// TEST(BulkSphere, DISABLED_HandlesInput) {
+TEST(BulkSphere, MultiLayerDu) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
   for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
@@ -81,7 +81,7 @@ TEST(BulkSphere, DISABLED_HandlesInput) {
 }
 
 //******************************************************************************
-TEST(BulkSphere, MesoMie) {
+TEST(BulkSphere, MesoMieDu) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
   nmie::MesoMie<nmie::FloatType> mesomie;
   for (const auto& data : parameters_and_results) {
@@ -94,7 +94,6 @@ TEST(BulkSphere, MesoMie) {
                     m * m,    // eps_m
                     {0, 0},   // d_parallel
                     {0, 0});  // d_perp
-                              // eps_m * xd / (eps_d * xm)
     mesomie.calc_Q();
 
     double Qext = static_cast<double>(mesomie.GetQext());

tests/test_py.py

@@ -1,6 +1,6 @@
 import unittest
-import numpy as np
-from scattnlay import scattnlay, mie, mie_mp
+from scattnlay import mie, mie_mp
+from scattnlay import mesomie
 
 test_cases = [
     # x, {Re(m), Im(m)}, Qext, Qsca, test_name
@@ -22,11 +22,36 @@ test_cases = [
 
 class TestStringMethods(unittest.TestCase):
 
-    def test_bulk(self):
+    def test_bulk_mesomie(self):
+        tol = 3e-7
+        for solver in [mesomie]:
+            if solver is None:
+                continue
+            print('Using solver: ', solver)
+            for case in test_cases:
+                x = case[0]
+                m = case[1]
+                solver.calc_ab(x,        # R
+                               x,   # xd
+                               x * m,    # xm
+                               1,   # eps_d
+                               m * m,    # eps_m
+                               0,   # d_parallel
+                               0)   # d_perp
+                # solver.calc_ab()
+                solver.calc_Q()
+                Qext = solver.GetQext()
+                Qsca = solver.GetQsca()
+                # print(Qext)
+                self.assertTrue((case[2]-Qext)/Qext < tol)
+                self.assertTrue((case[3]-Qsca)/Qsca < tol)
+
+    def test_bulk_multilayer(self):
         tol = 3e-7
         for solver in [mie, mie_mp]:
             if solver is None:
                 continue
+            print('Using solver: ', solver)
             for case in test_cases:
                 solver.SetLayersSize(case[0])
                 solver.SetLayersIndex(case[1])