get changes from component branch to compile standalone

src/nearfield.cc

@@ -241,7 +241,7 @@ int main(int argc, char *argv[]) {
       }
     }
 
-    nmie::nField(L, -1, x, m, -1, total_points, Xp, Yp, Zp, E, H);
+    nmie::nField(L, -1, x, m, -1, nmie::Modes::kAll, nmie::Modes::kAll, total_points, Xp, Yp, Zp, E, H);
 
     if (has_comment)
       printf("%6s\n", comment.c_str());

src/nmie-impl.hpp

@@ -832,6 +832,7 @@ namespace nmie {
     //      http://en.wikipedia.org/wiki/Loss_of_significance
     for (int i = nmax_ - 2; i >= 0; i--) {
       const int n = i + 1;
+      if (mode_n_ != Modes::kAll && n != mode_n_) continue;
       // Equation (27)
       Qext_ += (n + n + 1.0)*(an_[i].real() + bn_[i].real());
       // Equation (28)
@@ -1007,16 +1008,28 @@ namespace nmie {
   //   Rho: Radial distance                                                           //
   //   Phi: Azimuthal angle                                                           //
   //   Theta: Polar angle                                                             //
+  //   mode_n: mode order.                                                            //
+  //          -1 - use all modes (all_)                                               //
+  //           1 - use dipole mode only                                               //
+  //           2 - use quadrupole mode only                                           //
+  //           ...                                                                    //
+  //   mode_type: only used when mode_n != -1                                         //
+  //          0 - electric only                                                       //
+  //          1 - magnetic only                                                       //
+  //                                                                                  //
   //                                                                                  //
   // Output parameters:                                                               //
   //   E, H: Complex electric and magnetic fields                                     //
   //**********************************************************************************//
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcField(const FloatType Rho, const FloatType Theta, const FloatType Phi,
-                                std::vector<std::complex<FloatType> >& E, std::vector<std::complex<FloatType> >& H)  {
+  void MultiLayerMie<FloatType>::calcFieldByComponents(const FloatType Rho,
+                                const FloatType Theta, const FloatType Phi,
+                                std::vector<std::complex<FloatType> >& E,
+                                std::vector<std::complex<FloatType> >& H)  {
 
     std::complex<FloatType> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
-    std::vector<std::complex<FloatType> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
+    // Vector containing precomputed integer powers of i to avoid computation
+    std::vector<std::complex<FloatType> > ipow = {c_one, c_i, -c_one, -c_i};
     std::vector<std::complex<FloatType> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
     std::vector<std::complex<FloatType> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
     std::vector<std::complex<FloatType> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
@@ -1063,12 +1076,35 @@ namespace nmie {
       std::complex<FloatType> En = ipow[n1 % 4]
       *static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
       for (int i = 0; i < 3; i++) {
-        // electric field E [V m - 1] = EF*E0
-        E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
-              + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
-
-        H[i] += En*(-dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
-              +  c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
+        if (mode_n_ == Modes::kAll) {
+          // electric field E [V m - 1] = EF*E0
+          E[i] += En*(      cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
+                      + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
+
+          H[i] += En*(     -dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
+                      + c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
+          continue;
+        }
+        if (n1 == mode_n_) {
+          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]);
+
+            H[i] += En*(-dln_[l][n]*M1e1n[i]
+                        +aln_[l][n]*M3e1n[i]);
+            //std::cout << mode_n_;
+          }
+          if (mode_type_ == Modes::kMagnetic  || mode_type_ == Modes::kAll) {
+            E[i] += En*(  cln_[l][n]*M1o1n[i]
+                        - bln_[l][n]*M3o1n[i]);
+
+            H[i] += En*( -c_i*cln_[l][n]*N1o1n[i]
+                        + c_i*bln_[l][n]*N3o1n[i]);
+            //std::cout << mode_n_;
+          }
+          //std::cout << std::endl;
+        }
+        //throw std::invalid_argument("Error! Unexpected mode for field evaluation!\n mode_n="+std::to_string(mode_n)+", mode_type="+std::to_string(mode_type)+"\n=====*****=====");
       }
     }  // end of for all n
 
@@ -1077,7 +1113,7 @@ namespace nmie {
     for (int i = 0; i < 3; i++) {
       H[i] = hffact*H[i];
     }
-   }  // end of MultiLayerMie::calcField(...)
+   }  // end of MultiLayerMie::calcFieldByComponents(...)
 
 
   //**********************************************************************************//
@@ -1095,6 +1131,14 @@ namespace nmie {
   //   ncoord: Number of coordinate points                                            //
   //   Coords: Array containing all coordinates where the complex electric and        //
   //           magnetic fields will be calculated                                     //
+  //   mode_n: mode order.                                                            //
+  //          -1 - use all modes (all_)                                               //
+  //           1 - use dipole mode only                                               //
+  //           2 - use quadrupole mode only                                           //
+  //           ...                                                                    //
+  //   mode_type: only used when mode_n != -1                                         //
+  //          0 - electric only                                                       //
+  //          1 - magnetic only                                                       //
   //                                                                                  //
   // Output parameters:                                                               //
   //   E, H: Complex electric and magnetic field at the provided coordinates          //
@@ -1133,12 +1177,9 @@ namespace nmie {
       // If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
       Theta = (Rho > 0.0) ? nmm::acos(Zp/Rho) : 0.0;
 
-      // If Xp=Yp=0 then Phi is undefined. Just set it to zero to avoid problems
-      if (Xp == 0.0)
-        Phi = (Yp != 0.0) ? nmm::asin(Yp/nmm::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
-      else
-        Phi = nmm::acos(Xp/nmm::sqrt(pow2(Xp) + pow2(Yp)));
-      if (Yp < 0.0) Phi *= -1;
+      // std::atan2 should take care of any special cases, e.g.  Xp=Yp=0, etc.
+      Phi = nmm::atan2(Yp,Xp);
+
       // Avoid convergence problems due to Rho too small
       if (Rho < 1e-5) Rho = 1e-5;
       // std::cout << "Xp: "<<Xp<< "  Yp: "<<Yp<< "  Zp: "<<Zp<<std::endl;
@@ -1154,7 +1195,7 @@ namespace nmie {
       std::vector<std::complex<FloatType> > Es(3), Hs(3);
 
       // Do the actual calculation of electric and magnetic field
-      calcField(Rho, Theta, Phi, Es, Hs);
+      calcFieldByComponents(Rho, Theta, Phi, Es, Hs);
       for (int sph_coord = 0; sph_coord<3; ++sph_coord) {
         Es_[point][sph_coord] = Es[sph_coord];
         Hs_[point][sph_coord] = Hs[sph_coord];

src/nmie.cc

@@ -106,6 +106,54 @@ namespace nmie {
   }
 
   //**********************************************************************************//
+  // This function emulates a C call to calculate the scattering 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]     //
+  //   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:                                                               //
+  //   aln, bln ,cln, dln : Complex scattering amplitudes                             //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations in each layer   //
+  //**********************************************************************************//
+  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> > >& aln, std::vector<std::vector<std::complex<double> > >& bln, std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
+
+    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);
+
+      ml_mie.calcScattCoeffs();
+      ml_mie.calcExpanCoeffs();
+
+      aln = ConvertComplexVectorVector<double>(ml_mie.GetAln());
+      bln = ConvertComplexVectorVector<double>(ml_mie.GetBln());
+      cln = ConvertComplexVectorVector<double>(ml_mie.GetCln());
+      dln = ConvertComplexVectorVector<double>(ml_mie.GetDln());
+
+      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    //
   // and amplitudes.                                                                  //
   //                                                                                  //
@@ -137,7 +185,8 @@ namespace nmie {
   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) {
+           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!");
@@ -150,6 +199,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();
 
@@ -184,6 +234,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:                                                                //
@@ -309,8 +394,9 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  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 unsigned int ncoord,
+  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_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec,
              std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
     if (x.size() != L || m.size() != L)
@@ -332,6 +418,9 @@ namespace nmie {
       ml_mie.SetFieldCoords({ConvertVector<FloatType>(Xp_vec),
 	    ConvertVector<FloatType>(Yp_vec),
 	    ConvertVector<FloatType>(Zp_vec) });
+
+      ml_mie.SetModeNmaxAndType(mode_n, mode_type);
+
       ml_mie.RunFieldCalculation();
       E = ConvertComplexVectorVector<double>(ml_mie.GetFieldE());
       H = ConvertComplexVectorVector<double>(ml_mie.GetFieldH());

src/nmie.hpp

@@ -31,7 +31,7 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#define VERSION "2.2"
+#define VERSION "2.2"  //Compare with Makefile and setup.py
 #include <array>
 #include <complex>
 #include <cstdlib>
@@ -42,11 +42,58 @@
 #endif
 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 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);
-  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 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);
+  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);
+  // 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 per mode evaluation
+  enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};
 
   template <typename FloatType = double>
   class MultiLayerMie {    
@@ -61,6 +108,7 @@ namespace nmie {
     void RunMieCalculation();
     void RunFieldCalculation();
     void calcScattCoeffs();
+    void calcExpanCoeffs();
 
     // Return calculation results
     FloatType GetQext();
@@ -76,6 +124,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> > > GetAln(){return aln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetBln(){return bln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetCln(){return cln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetDln(){return dln_;};
+
     // Problem definition
     // Modify size of all layers
     void SetLayersSize(const std::vector<FloatType>& layer_size);
@@ -87,6 +140,8 @@ 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);
@@ -126,7 +181,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_;
 
@@ -158,8 +212,9 @@ namespace nmie {
                        std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n, 
                        std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n);
 
-    void calcField(const FloatType Rho, const FloatType Theta, const FloatType Phi,
-                   std::vector<std::complex<FloatType> >& E, std::vector<std::complex<FloatType> >& H);
+    void calcFieldByComponents(const FloatType Rho, const FloatType Theta, const FloatType Phi,
+                               std::vector<std::complex<FloatType> >& E,
+                               std::vector<std::complex<FloatType> >& H);
 
     bool isExpCoeffsCalc_ = false;
     bool isScaCoeffsCalc_ = false;
@@ -170,6 +225,9 @@ namespace nmie {
     // 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;
     int nmax_preset_ = -1;