Improved function to return expansion coefficients.

scattnlay/__init__.py

@@ -30,4 +30,4 @@
 #    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, scattnlay, fieldnlay
+from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay

scattnlay/main.py

@@ -126,7 +126,11 @@ def expancoeffs(x, m, li=-1, 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 expancoeffs_(x, m, nmax=nmax, pl=pl)
+            terms, an, bn, cn, dn = expancoeffs_(x, m, nmax=nmax, pl=pl)
+            if li >= 0 and li <= x.shape:
+                return terms, an[li], bn[li], cn[li], dn[li]
+            else:
+                return terms, an, bn, cn, dn
         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:
@@ -142,27 +146,30 @@ def expancoeffs(x, m, li=-1, nmax=-1, pl=-1, mp=False):
         nstore = nmax
 
     terms = np.zeros((x.shape[0]), dtype=int)
-    an = np.zeros((0, nstore), dtype=complex)
-    bn = np.zeros((0, nstore), dtype=complex)
-    cn = np.zeros((0, nstore), dtype=complex)
-    dn = np.zeros((0, nstore), dtype=complex)
+    an = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    bn = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    cn = np.zeros((0, x.shape[1], nstore), dtype=complex)
+    dn = np.zeros((0, x.shape[1], nstore), dtype=complex)
 
     for i, xi in enumerate(x):
         terms[i], a, b, c, d = expancoeffs_(xi, m[i], li=li, nmax=nmax, pl=pl)
 
         if terms[i] > nstore:
             nstore = terms[i]
-            an.resize((an.shape[0], nstore))
-            bn.resize((bn.shape[0], nstore))
-            cn.resize((cn.shape[0], nstore))
-            dn.resize((dn.shape[0], nstore))
+            an.resize((an.shape[0], an.shape[1], nstore))
+            bn.resize((bn.shape[0], bn.shape[1], nstore))
+            cn.resize((cn.shape[0], cn.shape[1], nstore))
+            dn.resize((dn.shape[0], dn.shape[1], nstore))
 
         an = np.vstack((an, a))
         bn = np.vstack((bn, b))
         cn = np.vstack((cn, c))
         dn = np.vstack((dn, d))
 
-    return terms, an, bn, cn, dn
+    if li >= 0 and li <= x.shape:
+        return terms, an[:, li, :], bn[:, li, :], cn[:, li, :], dn[:, li, :]
+    else:
+        return terms, an, bn, cn, dn
 #expancoeffs()
 
 

src/nmie-pybind11.cc

@@ -72,7 +72,7 @@ std::vector<T> flatten(const std::vector<std::vector<T>>& v) {
 
 
 template <typename T>
-py::array VectorVector2Py(const std::vector<std::vector<T > > &x) {
+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);
@@ -119,17 +119,19 @@ py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::arra
 
 py::tuple py_ExpanCoeffs(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 li=-1, const int nmax=-1, const int pl=-1) {
+                         const int nmax=-1, const int pl=-1) {
 
   auto c_x = Py2Vector<double>(py_x);
   auto c_m = Py2Vector< std::complex<double> >(py_m);
 
   int terms = 0;
-  std::vector<std::complex<double> > c_an, c_bn, c_cn, c_dn;
   int L = py_x.size();
-  terms = nmie::ExpanCoeffs(L, pl, c_x, c_m, li, nmax, c_an, c_bn, c_cn, c_dn);
+
+  std::vector<std::vector<std::complex<double> > > c_an(L + 1), c_bn(L + 1), c_cn(L + 1), c_dn(L + 1);
+
+  terms = nmie::ExpanCoeffs(L, pl, c_x, c_m, nmax, c_an, c_bn, c_cn, c_dn);
   
-  return py::make_tuple(terms, VectorComplex2Py(c_an), VectorComplex2Py(c_bn), VectorComplex2Py(c_cn), VectorComplex2Py(c_dn));
+  return py::make_tuple(terms, Vector2DComplex2Py(c_an), Vector2DComplex2Py(c_bn), Vector2DComplex2Py(c_cn), Vector2DComplex2Py(c_dn));
 }
 
 
@@ -172,8 +174,8 @@ py::tuple py_fieldnlay(const py::array_t<double, py::array::c_style | py::array:
   for (auto& f : H) f.resize(3);
   int L = py_x.size(), terms;
   terms = nmie::nField(L, pl, c_x, c_m, nmax, ncoord, c_Xp, c_Yp, c_Zp, E, H);
-  auto py_E = VectorVector2Py<std::complex<double> >(E);
-  auto py_H = VectorVector2Py<std::complex<double> >(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-pybind11.hpp

@@ -45,7 +45,7 @@ py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::arra
 
 py::tuple py_ExpanCoeffs(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 li=-1, const int nmax=-1, const int pl=-1);
+                         const int nmax=-1, const int pl=-1);
 
 
 py::tuple py_scattnlay(const py::array_t<double, py::array::c_style | py::array::forcecast> &py_x,

src/nmie.cc

@@ -104,7 +104,7 @@ namespace nmie {
     }
     return 0;
   }
-  
+
   //**********************************************************************************//
   // This function emulates a C call to calculate the expansion coefficients          //
   // required to calculate both the near- and far-field parameters.                   //
@@ -114,21 +114,19 @@ namespace nmie {
   //   pl: Index of PEC layer. If there is none just send -1                          //
   //   x: Array containing the size parameters of the layers [0..L-1]                 //
   //   m: Array containing the relative refractive indexes of the layers [0..L-1]     //
-  //   li: Index of layer to get expansion coefficients. -1 means outside the sphere  //
   //   nmax: Maximum number of multipolar expansion terms to be used for the          //
   //         calculations. Only use it if you know what you are doing, otherwise      //
   //         set this parameter to -1 and the function will calculate it.             //
   //                                                                                  //
   // Output parameters:                                                               //
-  //   an[li], bn[li], cn[li], dn[li]: Complex expansion coefficients of ith layer    //
+  //   an, bn, cn, dn: Complex expansion coefficients                                 //
   //                                                                                  //
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
-                  const int li, const int nmax,
-                  std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn,
-                  std::vector<std::complex<double> >& cn, std::vector<std::complex<double> >& dn) {
+  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+                  std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn,
+                  std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -144,17 +142,10 @@ namespace nmie {
     // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
       ml_mie.calcExpanCoeffs();
 
-      if (li >= L || li < 0) { // Just return the scattering coefficients
-        an = ConvertComplexVector<double>(ml_mie.GetLayerAn()[L]);
-        bn = ConvertComplexVector<double>(ml_mie.GetLayerBn()[L]);
-        cn = ConvertComplexVector<double>(ml_mie.GetLayerCn()[L]);
-        dn = ConvertComplexVector<double>(ml_mie.GetLayerDn()[L]);
-      } else { // Return the expansion coefficients for ith layer
-        an = ConvertComplexVector<double>(ml_mie.GetLayerAn()[li]);
-        bn = ConvertComplexVector<double>(ml_mie.GetLayerBn()[li]);
-        cn = ConvertComplexVector<double>(ml_mie.GetLayerCn()[li]);
-        dn = ConvertComplexVector<double>(ml_mie.GetLayerDn()[li]);
-      }
+      an = ConvertComplexVectorVector<double>(ml_mie.GetLayerAn());
+      bn = ConvertComplexVectorVector<double>(ml_mie.GetLayerBn());
+      cn = ConvertComplexVectorVector<double>(ml_mie.GetLayerCn());
+      dn = ConvertComplexVectorVector<double>(ml_mie.GetLayerDn());
 
       return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {

src/nmie.hpp

@@ -40,7 +40,7 @@
 #include <boost/math/constants/constants.hpp>
 namespace nmie {
   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 ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, const int li, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn, std::vector<std::complex<double> >& cn, std::vector<std::complex<double> >& dn);
+  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn, std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn);
   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, 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);