adding mode selection into scattnlay

examples/fieldplot.py

@@ -277,17 +277,17 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
                         )
         ax.axis("image")
 
-        # Add colorbar
-        cbar = fig.colorbar(cax, ticks=[a for a in scale_ticks], ax=ax
-                            #,fraction=0.45
-            )
-        # vertically oriented colorbar
-        if 'angle' in field_to_plot:
-            cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
-        else:
-            cbar.ax.set_yticklabels(['%g' % (a) for a in scale_ticks])
-        # pos = list(cbar.ax.get_position().bounds)
-        #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
+        # # Add colorbar
+        # cbar = fig.colorbar(cax, ticks=[a for a in scale_ticks], ax=ax
+        #                     #,fraction=0.45
+        #     )
+        # # vertically oriented colorbar
+        # if 'angle' in field_to_plot:
+        #     cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
+        # else:
+        #     cbar.ax.set_yticklabels(['%g' % (a) for a in scale_ticks])
+        # # pos = list(cbar.ax.get_position().bounds)
+        # #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
         lp2 = -10.0
         lp1 = -1.0
         if crossplane == 'XZ':

examples/force-quad.py

@@ -50,7 +50,7 @@ WL = 550
 x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
 m[0] = index_Ag/nm
 
-R_st = x[0]*4.11
+R_st = x[0]*2
 #R_st = 0.31
 
 dx = R_st*4.0
@@ -100,7 +100,7 @@ def dipole(coord):
 
 def force(in_coord):
     coord = in_coord.T
-    terms, Eall, Hall = fieldnlay(np.array([x]), np.array([m]), coord, mode_n=-1, mode_type=0)
+    terms, Eall, Hall = fieldnlay(np.array([x]), np.array([m]), coord)#, mode_n=-1, mode_type=0)
     E_all = Eall[0, :, :]
     H_all = Hall[0, :, :]
 

scattnlay.pyx

@@ -54,7 +54,7 @@ cdef extern from "py_nmie.h":
     #                 vector[vector[double] ]&  dlnr,
     #                 vector[vector[double] ]&  dlni)
     cdef int ExpansionCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double alnr[], double alni[], double blnr[], double blni[], double clnr[], double clni[], double dlnr[], double dlni[])
-    cdef int nMie(int L, int pl, vector[double] x, vector[complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[])
+    cdef int nMie(int L, int pl, vector[double] x, vector[complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[], int mode_n, int mode_type)
     cdef int nField(int L, int pl, vector[double] x, vector[complex] m, int nmax, int mode_n, int mode_type, int nCoords, vector[double] Xp, vector[double] Yp, vector[double] Zp, double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[], double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[])
 
 def scattcoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.int_t nmax, np.int_t pl = -1):
@@ -123,7 +123,7 @@ def expansioncoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex1
             dln[i][l] = dlnr[l*nmax:(l+1)*nmax].copy('C') + 1.0j*dlni[l*nmax:(l+1)*nmax].copy('C')
     return terms, aln, bln, cln, dln
 
-def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t nmax = -1, np.int_t pl = -1):
+def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t nmax = -1, np.int_t mode_n = -1, np.int_t mode_type = -1, np.int_t pl = -1):
     cdef Py_ssize_t i
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
@@ -150,7 +150,7 @@ def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t,
         S2r = np.zeros(theta.shape[0], dtype = np.float64)
         S2i = np.zeros(theta.shape[0], dtype = np.float64)
 
-        terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i))
+        terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i), mode_n, mode_type)
 
         S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C')
         S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C')

src/nmie-impl.hpp

@@ -1085,7 +1085,7 @@ namespace nmie {
           continue;
         }
         if (n1 == mode_n) {
-          if (mode_type == Modes::kElectric) {
+          if (mode_type == Modes::kElectric || mode_type == Modes::kAll) {
             E[i] += En*( -c_i*dln_[l][n]*N1e1n[i]
                         + c_i*aln_[l][n]*N3e1n[i]);
 
@@ -1093,7 +1093,7 @@ namespace nmie {
                         +aln_[l][n]*M3e1n[i]);
             continue;
           }
-          if (mode_type == Modes::kMagnetic) {
+          if (mode_type == Modes::kMagnetic  || mode_type == Modes::kAll) {
             E[i] += En*(  cln_[l][n]*M1o1n[i]
                         - bln_[l][n]*M3o1n[i]);
 

src/nmie.cc

@@ -173,7 +173,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  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, 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 mode_n, int mode_type) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -186,6 +186,7 @@ namespace nmie {
       ml_mie.SetAngles(ConvertVector<FloatType>(Theta));
       ml_mie.SetPECLayer(pl);
       ml_mie.SetMaxTerms(nmax);
+      ml_mie.SetModeNmaxAndType(mode_n, mode_type);
 
       ml_mie.RunMieCalculation();
 
@@ -220,6 +221,41 @@ namespace nmie {
   // This function is just a wrapper to call the full 'nMie' function with fewer      //
   // parameters, it is here mainly for compatibility with older versions of the       //
   // program. Also, you can use it if you neither have a PEC layer nor want to define //
+  // any limit for the maximum number of terms nor limit to some mode.                //
+  //                                                                                  //
+  // 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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //   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:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  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) {
+    return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, -1, -1);
+
+  }
+  //**********************************************************************************//
+  // This function is just a wrapper to call the full 'nMie' function with fewer      //
+  // parameters, it is here mainly for compatibility with older versions of the       //
+  // program. Also, you can use it if you neither have a PEC layer nor want to define //
   // any limit for the maximum number of terms.                                       //
   //                                                                                  //
   // Input parameters:                                                                //

src/nmie.hpp

@@ -37,13 +37,51 @@
 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 ExpansionCoeffs(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);
-  int nMie(const unsigned int L, 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);
+  // pl, nmax, mode_n, mode_type
+  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 mode_n, int mode_type);
+  // pl and nmax
+  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);
+  // no pl and nmax
+  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);
+  // pl
+  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);
+  // nmax
+  int nMie(const unsigned int L,
+           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 nField(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const int mode_n, const int mode_type, const unsigned int ncoord, const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
 
-  // constants for field evaluation
+  // constants for per mode evaluation
   enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};
 
   template <typename FloatType = double>
@@ -91,7 +129,9 @@ namespace nmie {
     void SetFieldCoords(const std::vector< std::vector<FloatType> >& coords);
     // Modify index of PEC layer
     void SetPECLayer(int layer_position = 0);
-
+    // Modify the mode taking into account for evaluation of output variables
+    void SetModeNmaxAndType(int mode_n, int mode_type){mode_n_ = mode_n; mode_type_ = mode_type;};
+    
     // Set a fixed value for the maximun number of terms
     void SetMaxTerms(int nmax);
     // Get maximun number of terms
@@ -173,6 +213,9 @@ namespace nmie {
     std::vector<FloatType> theta_;
     // Should be -1 if there is no PEC.
     int PEC_layer_position_ = -1;
+    
+    int mode_n_ = Modes::kAll;
+    int mode_type_ = Modes::kAll;
 
     // with calcNmax(int first_layer);
     int nmax_ = -1;

src/py_nmie.cc

@@ -111,17 +111,20 @@ int ExpansionCoeffs(const unsigned int L, const int pl,
 // Same as nMie in 'nmie.h' but uses double arrays to return the results (useful for python).
 // This is a workaround because I have not been able to return the results using 
 // std::vector<std::complex<double> >
-int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
+int nMie(const int L, const int pl,
+         std::vector<double>& x, std::vector<std::complex<double> >& m,
          const int nTheta, std::vector<double>& Theta, const int nmax,
-         double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-		 double S1r[], double S1i[], double S2r[], double S2i[]) {
+         double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
+         double *g, double *Albedo,
+         double S1r[], double S1i[], double S2r[], double S2i[],
+         int mode_n, int mode_type) {
 
   int i, result;
   std::vector<std::complex<double> > S1, S2;
   S1.resize(nTheta);
   S2.resize(nTheta);
 
-  result = nmie::nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+  result = nmie::nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, mode_n, mode_type);
 
   for (i = 0; i < nTheta; i++) {
     S1r[i] = S1[i].real();

src/py_nmie.h

@@ -51,10 +51,13 @@ int ExpansionCoeffs(const unsigned int L, const int pl,
                     double alnr[], double alni[], double blnr[], double blni[],
                     double clnr[], double clni[], double dlnr[], double dlni[]);
 
-int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
+int nMie(const int L, const int pl,
+         std::vector<double>& x, std::vector<std::complex<double> >& m,
          const int nTheta, std::vector<double>& Theta, const int nmax,
-         double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-		 double S1r[], double S1i[], double S2r[], double S2i[]);
+         double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
+         double *g, double *Albedo,
+         double S1r[], double S1i[], double S2r[], double S2i[],
+         int mode_n, int mode_type);
 
 int nField(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, const int mode_n, const int mode_type,
            const int nCoords, std::vector<double>& Xp, std::vector<double>& Yp, std::vector<double>& Zp,