Implemented Python function to return expansion coefficients.

scattnlay/main.py

@@ -68,6 +68,10 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
     elif len(x.shape) != 2:
         raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     if nmax == -1:
         nstore = 0
     else:
@@ -78,23 +82,88 @@ 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)
+
+        if terms[i] > nstore:
+            nstore = terms[i]
+            an.resize((an.shape[0], nstore))
+            bn.resize((bn.shape[0], nstore))
+
+        an = np.vstack((an, a))
+        bn = np.vstack((bn, b))
+
+    return terms, an, bn
+#scattcoeffs()
+
+def expancoeffs(x, m, li=-1, nmax=-1, pl=-1, mp=False):
+    """
+    expancoeffs(x, m[, li, nmax, pl, mp])
+
+    Calculate the scattering coefficients required to calculate both the
+    near- and far-field parameters.
+
+        x: Size parameters (1D or 2D ndarray)
+        m: Relative refractive indices (1D or 2D ndarray)
+        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.
+        pl: Index of PEC layer. If there is none just send -1.
+        mp: Use multiple (True) or double (False) precision.
+
+    Returns: (terms, an, bn, cn, dn)
+    with
+        terms: Number of multipolar expansion terms used for the calculations
+        an, bn, cn, dn: Complex expansion coefficients of ith layer
+    """
+
+    if mp:
+        from scattnlay_mp import expancoeffs as expancoeffs_
+    else:
+        from scattnlay_dp import expancoeffs as expancoeffs_
+
+    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:
-            mi = m
+            return expancoeffs_(x, m, nmax=nmax, pl=pl)
         else:
-            mi = m[i]
+            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:
+        raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
+
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
+    if nmax == -1:
+        nstore = 0
+    else:
+        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)
 
-        terms[i], a, b = scattcoeffs_(xi, mi, nmax=nmax, pl=pl)
+    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 = np.vstack((an, a))
         bn = np.vstack((bn, b))
+        cn = np.vstack((cn, c))
+        dn = np.vstack((dn, d))
 
-    return terms, an, bn
-#scattcoeffs()
+    return terms, an, bn, cn, dn
+#expancoeffs()
 
 
 def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
@@ -142,6 +211,10 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     if len(theta.shape) != 1:
         raise ValueError('The scattering angles (theta) should be a 1-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     terms = np.zeros((x.shape[0]), dtype=int)
     Qext = np.zeros((x.shape[0]), dtype=float)
     Qsca = np.zeros((x.shape[0]), dtype=float)
@@ -154,12 +227,7 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
     S2 = np.zeros((x.shape[0], theta.shape[0]), dtype=complex)
 
     for i, xi in enumerate(x):
-        if len(m.shape) == 1:
-            mi = m
-        else:
-            mi = m[i]
-
-        terms[i], Qext[i], Qsca[i], Qabs[i], Qbk[i], Qpr[i], g[i], Albedo[i], S1[i], S2[i] = scattnlay_(xi, mi, theta, nmax=nmax, pl=pl)
+        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
 #scattnlay()
@@ -174,11 +242,11 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
         x: Size parameters (1D or 2D ndarray)
         m: Relative refractive indices (1D or 2D ndarray)
         xp: Array containing all X coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         yp: Array containing all Y coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         zp: Array containing all Z coordinates to calculate the complex
-            electric and magnetic fields (1D ndarray)
+            electric and magnetic fields (1D* ndarray)
         nmax: Maximum number of multipolar expansion terms to be used for the
               calculations. Only use it if you know what you are doing.
         pl: Index of PEC layer. If there is none just send -1.
@@ -188,6 +256,9 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     with
         terms: Number of multipolar expansion terms used for the calculations
         E, H: Complex electric and magnetic field at the provided coordinates
+
+    *Note: We assume that the coordinates are referred to the first wavelength
+           (or structure) and correct it for the following ones
     """
 
     if mp:
@@ -205,17 +276,18 @@ def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
     elif len(x.shape) != 2:
         raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
 
+    # Repeat the same m for all wavelengths
+    if len(m.shape) == 1:
+        m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
+
     terms = np.zeros((x.shape[0]), dtype=int)
     E = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
     H = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
 
     for i, xi in enumerate(x):
-        if len(m.shape) == 1:
-            mi = m
-        else:
-            mi = m[i]
-
-        terms[i], E[i], H[i] = fieldnlay_(xi, mi, xp, yp, zp, nmax=nmax, pl=pl)
+        # (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)
 
     return terms, E, H
 #fieldnlay()

src/nmie-pybind11.cc

@@ -117,6 +117,22 @@ 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) {
+
+  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);
+  
+  return py::make_tuple(terms, VectorComplex2Py(c_an), VectorComplex2Py(c_bn), VectorComplex2Py(c_cn), VectorComplex2Py(c_dn));
+}
+
+
 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,

src/nmie-pybind11.hpp

@@ -43,6 +43,11 @@ py::tuple py_ScattCoeffs(const py::array_t<double, py::array::c_style | py::arra
                          const int nmax=-1, const int pl=-1);
 
 
+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);
+
+
 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,

src/nmie.cc

@@ -104,6 +104,67 @@ 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.                   //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   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    //
+  //                                                                                  //
+  // 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) {
+
+    if (x.size() != L || m.size() != L)
+        throw std::invalid_argument("Declared number of layers do not fit x and m!");
+    try {
+      MultiLayerMie<FloatType> ml_mie;
+      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
+      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
+      ml_mie.SetPECLayer(pl);
+      ml_mie.SetMaxTerms(nmax);
+
+    // Calculate scattering coefficients an_ and bn_
+      ml_mie.calcScattCoeffs();
+    // 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]);
+      }
+
+      return ml_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;
+      throw std::invalid_argument(ia);
+      return -1;
+    }
+    return 0;
+  }
 
   //**********************************************************************************//
   // This function emulates a C call to calculate the actual scattering parameters    //

src/nmie.hpp

@@ -40,6 +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 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);
@@ -59,6 +60,7 @@ namespace nmie {
     void RunMieCalculation();
     void RunFieldCalculation();
     void calcScattCoeffs();
+    void calcExpanCoeffs();
 
     // Return calculation results
     FloatType GetQext();
@@ -74,6 +76,11 @@ namespace nmie {
     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_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerBn(){return bln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerCn(){return cln_;};
+    std::vector< std::vector<std::complex<FloatType> > > GetLayerDn(){return dln_;};
+
     // Problem definition
     // Modify size of all layers
     void SetLayersSize(const std::vector<FloatType>& layer_size);
@@ -124,7 +131,6 @@ namespace nmie {
     // Scattering coefficients
     std::vector<std::complex<FloatType> > an_, bn_;
     std::vector< std::vector<std::complex<FloatType> > > aln_, bln_, cln_, dln_;
-    void calcExpanCoeffs();
     // Points for field evaluation
     std::vector< std::vector<FloatType> > coords_;