update

examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -40,11 +40,11 @@ print("Using scattnlay from ", inspect.getfile(scattnlay))
 
 npts = 351
 # npts = 11
-factor = 3 # plot extent compared to sphere radius
+factor = 2 # plot extent compared to sphere radius
 index_H2O = 1.33+(1e-6)*1j
 
 WL = 0.532 #mkm
-total_r = 1 #mkm
+total_r = 125 #mkm
 isMP = False
 # isMP = True
 
@@ -89,14 +89,14 @@ Er = np.absolute(Ec)
 Eabs2 = (Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
 Eabs_data = np.resize(Eabs2, (npts, npts))
 label = r'$|E|^2$'
-pos = plt.imshow(Eabs_data,
+pos = plt.imshow(Eabs_data.T,
            cmap='gnuplot',
                  # cmap='jet',
            vmin=0., vmax=14
 
            )
 plt.colorbar(pos)
-print(np.min(Eabs_data), np.max(Eabs_data)," terms = "+str(terms))
+print(np.min(Eabs_data), np.max(Eabs_data)," terms = "+str(terms), ' size=', Eabs_data.size)
 mp = ''
 if isMP: mp = '_mp'
 plt.savefig("R"+str(total_r)+"mkm"+mp+".jpg",

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, expancoeffs, scattnlay, fieldnlay
+from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay, mie, mie_mp

scattnlay/main.py

@@ -31,11 +31,17 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 import numpy as np
 
+from scattnlay_mp import mie_mp as mie_mp_
+from scattnlay_dp import mie_dp
+
+mie = mie_dp()
+mie_mp = mie_mp_()
+
 def scattcoeffs_(x, m, nmax=-1, pl=-1, mp=False):
     if mp:
         from scattnlay_mp import mie_mp as mie_
     else:
-        from scattnlay3 import mie_dp as mie_
+        from scattnlay_dp import mie_dp as mie_
         # from scattnlay_mp import mie_mp as mie_
     mie = mie_()
     mie.SetLayersSize(x)

setup.py

@@ -65,10 +65,19 @@ For details see: O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354."""
                              ["src/pb11_wrapper.cc"],
                              language="c++",
                              include_dirs=[np.get_include(), pb.get_include()],
-                             extra_compile_args=['-std=c++11']),
+                             # extra_compile_args=['-std=c++11']),
+                             extra_compile_args=['-std=c++11', '-O3',
+                                                 '-mavx2', '-mfma',
+                                                 '-finline-limit=1000000', '-ffp-contract=fast']),
+
                    Extension("scattnlay_mp",
                              ["src/pb11_wrapper.cc"],
                              language="c++",
                              include_dirs=[np.get_include(), pb.get_include()],
-                             extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
+                             extra_compile_args=['-std=c++11', '-O3',
+                                                 '-mavx2', '-mfma',
+                                                 '-finline-limit=1000000', '-ffp-contract=fast',
+                                                 '-DMULTI_PRECISION=100']),
+                             # extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])
+                   ]
       )

src/nmie-basic.hpp

@@ -525,33 +525,36 @@ namespace nmie {
   // Output parameters:                                                               //
   //   Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics                     //
   //**********************************************************************************//
-  template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
-                                    const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
-                                    const FloatType &Pi, const FloatType &Tau, const FloatType &n,
-                                    std::vector<std::complex<FloatType> > &Mo1n, std::vector<std::complex<FloatType> > &Me1n,
-                                    std::vector<std::complex<FloatType> > &No1n, std::vector<std::complex<FloatType> > &Ne1n) {
+  template <typename FloatType>  template <typename evalType>
+  void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<evalType> Rho, const evalType Theta, const evalType Phi,
+                                    const std::complex<evalType> &rn, const std::complex<evalType> &Dn,
+                                    const evalType &Pi, const evalType &Tau, const evalType &n,
+                                    std::vector<std::complex<evalType> > &Mo1n, std::vector<std::complex<evalType> > &Me1n,
+                                    std::vector<std::complex<evalType> > &No1n, std::vector<std::complex<evalType> > &Ne1n) {
 
     // using eq 4.50 in BH
-    std::complex<FloatType> c_zero(0.0, 0.0);
+    std::complex<evalType> c_zero(0.0, 0.0);
 
     using nmm::sin;
     using nmm::cos;
+    auto sin_Phi = sin(Phi);
+    auto cos_Phi = cos(Phi);
+    auto sin_Theta = sin(Theta);
     Mo1n[0] = c_zero;
-    Mo1n[1] = cos(Phi)*Pi*rn/Rho;
-    Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
+    Mo1n[1] = cos_Phi*Pi*rn/Rho;
+    Mo1n[2] = -sin_Phi*Tau*rn/Rho;
 
     Me1n[0] = c_zero;
-    Me1n[1] = -sin(Phi)*Pi*rn/Rho;
-    Me1n[2] = -cos(Phi)*Tau*rn/Rho;
+    Me1n[1] = -sin_Phi*Pi*rn/Rho;
+    Me1n[2] = -cos_Phi*Tau*rn/Rho;
 
-    No1n[0] = sin(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
-    No1n[1] = sin(Phi)*Tau*Dn*rn/Rho;
-    No1n[2] = cos(Phi)*Pi*Dn*rn/Rho;
+    No1n[0] = sin_Phi*(n*n + n)*sin_Theta*Pi*rn/Rho/Rho;
+    No1n[1] = sin_Phi*Tau*Dn*rn/Rho;
+    No1n[2] = cos_Phi*Pi*Dn*rn/Rho;
 
-    Ne1n[0] = cos(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
-    Ne1n[1] = cos(Phi)*Tau*Dn*rn/Rho;
-    Ne1n[2] = -sin(Phi)*Pi*Dn*rn/Rho;
+    Ne1n[0] = cos_Phi*(n*n + n)*sin_Theta*Pi*rn/Rho/Rho;
+    Ne1n[1] = cos_Phi*Tau*Dn*rn/Rho;
+    Ne1n[2] = -sin_Phi*Pi*Dn*rn/Rho;
   }  // end of MultiLayerMie::calcSpherHarm(...)
 
 

src/nmie-nearfield.hpp

@@ -232,25 +232,25 @@ namespace nmie {
   // Output parameters:                                                               //
   //   E, H: Complex electric and magnetic fields                                     //
   //**********************************************************************************//
-  template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcFieldByComponents(const FloatType Rho,
-                                const FloatType Theta, const FloatType Phi,
-                                const std::vector<std::complex<FloatType> > &Psi,
-                                const std::vector<std::complex<FloatType> > &D1n,
-                                const std::vector<std::complex<FloatType> > &Zeta,
-                                const std::vector<std::complex<FloatType> > &D3n,
-                                const std::vector<FloatType> &Pi,
-                                const std::vector<FloatType> &Tau,
-                                std::vector<std::complex<FloatType> > &E,
-                                std::vector<std::complex<FloatType> > &H)  {
+  template <typename FloatType>  template <typename evalType>
+  void MultiLayerMie<FloatType>::calcFieldByComponents(const evalType Rho,
+                                const evalType Theta, const evalType Phi,
+                                const std::vector<std::complex<evalType> > &Psi,
+                                const std::vector<std::complex<evalType> > &D1n,
+                                const std::vector<std::complex<evalType> > &Zeta,
+                                const std::vector<std::complex<evalType> > &D3n,
+                                const std::vector<evalType> &Pi,
+                                const std::vector<evalType> &Tau,
+                                std::vector<std::complex<evalType> > &E,
+                                std::vector<std::complex<evalType> > &H)  {
     auto nmax = Psi.size() - 1;
-    std::complex<FloatType> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
+    std::complex<evalType> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
     // 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<evalType> > ipow = {c_one, c_i, -c_one, -c_i};
+    std::vector<std::complex<evalType> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
+    std::vector<std::complex<evalType> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
 
-    std::complex<FloatType> ml;
+    std::complex<evalType> ml;
 
     // Initialize E and H
     for (int i = 0; i < 3; i++) {
@@ -263,16 +263,16 @@ namespace nmie {
 
     for (unsigned int n = 0; n < nmax; n++) {
       int n1 = n + 1;
-      auto rn = static_cast<FloatType>(n1);
+      auto rn = static_cast<evalType>(n1);
 
       // using BH 4.12 and 4.50
       calcSpherHarm(Rho*ml, Theta, Phi, Psi[n1], D1n[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
       calcSpherHarm(Rho*ml, Theta, Phi, Zeta[n1], D3n[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
 
       // Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
-      std::complex<FloatType> En = ipow[n1 % 4]
-      *static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
-      std::complex<FloatType> Ediff, Hdiff;
+      std::complex<evalType> En = ipow[n1 % 4]
+      *static_cast<evalType>((rn + rn + 1.0)/(rn*rn + rn));
+      std::complex<evalType> Ediff, Hdiff;
       for (int i = 0; i < 3; i++) {
         Ediff = 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]);
@@ -318,7 +318,7 @@ namespace nmie {
     }  // end of for all n
 
     // magnetic field
-    std::complex<FloatType> hffact = ml/static_cast<FloatType>(cc_*mu_);
+    std::complex<evalType> hffact = ml/static_cast<evalType>(cc_*mu_);
     for (int i = 0; i < 3; i++) {
       H[i] = hffact*H[i];
     }
@@ -422,9 +422,9 @@ double eval_delta(const unsigned int steps, const double from_value, const doubl
 
 // ml - refractive index
 // l - Layer number
-template <typename FloatType>
-void MultiLayerMie<FloatType>::GetIndexAtRadius(const FloatType Rho,
-                                                std::complex<FloatType> &ml,
+template <typename FloatType> template <typename evalType>
+void MultiLayerMie<FloatType>::GetIndexAtRadius(const evalType Rho,
+                                                std::complex<evalType> &ml,
                                                 unsigned int &l) {
   l = 0;
   if (Rho > size_param_.back()) {
@@ -439,9 +439,9 @@ void MultiLayerMie<FloatType>::GetIndexAtRadius(const FloatType Rho,
     ml = refractive_index_[l];
   }
 }
-template <typename FloatType>
-void MultiLayerMie<FloatType>::GetIndexAtRadius(const FloatType Rho,
-                                                std::complex<FloatType> &ml) {
+template <typename FloatType> template <typename evalType>
+void MultiLayerMie<FloatType>::GetIndexAtRadius(const evalType Rho,
+                                                std::complex<evalType> &ml) {
   unsigned int l;
   GetIndexAtRadius(Rho, ml, l);
 }
@@ -506,10 +506,19 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int outer_arc_poin
       || outer_arc_points < 1 || radius_points < 1
       || from_Rho < 0.)
     throw std::invalid_argument("Error! Invalid argument for RunFieldCalculationPolar() !");
-
+  int theta_points = 0, phi_points = 0;
+  if (to_Theta-from_Theta > to_Phi-from_Phi) {
+    theta_points = outer_arc_points;
+    phi_points =  static_cast<int>((to_Phi-from_Phi)/(to_Theta-from_Theta) * outer_arc_points);
+  } else {
+    phi_points = outer_arc_points;
+    theta_points =  static_cast<int>((to_Theta-from_Theta)/(to_Phi-from_Phi) * outer_arc_points);
+  }
+  if (theta_points == 0) theta_points = 1;
+  if (phi_points == 0) phi_points = 1;
   calcMieSeriesNeededToConverge(to_Rho);
 
-  std::vector<std::vector<FloatType> >  Pi(outer_arc_points), Tau(outer_arc_points);
+  std::vector<std::vector<FloatType> >  Pi(theta_points), Tau(theta_points);
   calcPiTauAllTheta(from_Theta, to_Theta, Pi, Tau);
 
   std::vector<std::vector<std::complex<FloatType> > > Psi(radius_points), D1n(radius_points),
@@ -518,8 +527,8 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int outer_arc_poin
                                    Psi, D1n, Zeta, D3n);
 
   double delta_Rho = eval_delta<double>(radius_points, from_Rho, to_Rho);
-  double delta_Theta = eval_delta<double>(outer_arc_points, from_Theta, to_Theta);
-  double delta_Phi = eval_delta<double>(radius_points, from_Phi, to_Phi);
+  double delta_Theta = eval_delta<double>(theta_points, from_Theta, to_Theta);
+  double delta_Phi = eval_delta<double>(phi_points, from_Phi, to_Phi);
   Es_.clear(); Hs_.clear(); coords_polar_.clear();
   for (int j=0; j < radius_points; j++) {
     auto Rho = static_cast<FloatType>(from_Rho + j * delta_Rho);

src/nmie.hpp

@@ -204,7 +204,7 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
                                   const double from_Rho=0, const double to_Rho=static_cast<double>(1.),
                                   const double from_Theta=0, const double to_Theta=static_cast<double>(3.14159265358979323),
                                   const double from_Phi=0, const double to_Phi=static_cast<double>(3.14159265358979323),
-                                  const bool isIgnoreAvailableNmax = true); // TODO change to false for production
+                                  const bool isIgnoreAvailableNmax = false);
 
     void calcScattCoeffs();
     void calcExpanCoeffs();
@@ -245,8 +245,9 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
     void SetLayersIndex(const std::vector< std::complex<FloatType> > &index);
     template <typename inputType>
     void SetLayersIndex(const py::array_t<std::complex<inputType>, py::array::c_style | py::array::forcecast> &py_index);
-    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml, unsigned int &l);
-    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml);
+
+    template <typename evalType=FloatType> void GetIndexAtRadius(const evalType Rho, std::complex<evalType> &ml, unsigned int &l);
+    template <typename evalType=FloatType> void GetIndexAtRadius(const evalType Rho, std::complex<evalType> &ml);
     // Modify scattering (theta) py_angles
     void SetAngles(const std::vector<FloatType> &py_angles);
     template <typename inputType>
@@ -328,21 +329,23 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
                      std::vector<std::complex<FloatType> > &Zeta);
     void calcPiTau(const FloatType &costheta,
                    std::vector<FloatType> &Pi, std::vector<FloatType> &Tau);
-    void calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
-                       const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
-                       const FloatType &Pi, const FloatType &Tau, const FloatType &n,
-                       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 calcFieldByComponents(const FloatType Rho, const FloatType Theta, const FloatType Phi,
-                               const std::vector<std::complex<FloatType> > &Psi,
-                               const std::vector<std::complex<FloatType> > &D1n,
-                               const std::vector<std::complex<FloatType> > &Zeta,
-                               const std::vector<std::complex<FloatType> > &D3n,
-                               const std::vector<FloatType> &Pi,
-                               const std::vector<FloatType> &Tau,
-                               std::vector<std::complex<FloatType> > &E,
-                               std::vector<std::complex<FloatType> > &H);
+    template <typename evalType=FloatType>
+    void calcSpherHarm(const std::complex<evalType> Rho, const evalType Theta, const evalType Phi,
+                       const std::complex<evalType> &rn, const std::complex<evalType> &Dn,
+                       const evalType &Pi, const evalType &Tau, const evalType &n,
+                       std::vector<std::complex<evalType> > &Mo1n, std::vector<std::complex<evalType> > &Me1n,
+                       std::vector<std::complex<evalType> > &No1n, std::vector<std::complex<evalType> > &Ne1n);
+
+    template <typename evalType=FloatType>
+    void calcFieldByComponents(const evalType Rho, const evalType Theta, const evalType Phi,
+                               const std::vector<std::complex<evalType> > &Psi,
+                               const std::vector<std::complex<evalType> > &D1n,
+                               const std::vector<std::complex<evalType> > &Zeta,
+                               const std::vector<std::complex<evalType> > &D3n,
+                               const std::vector<evalType> &Pi,
+                               const std::vector<evalType> &Tau,
+                               std::vector<std::complex<evalType> > &E,
+                               std::vector<std::complex<evalType> > &H);
 
     bool isExpCoeffsCalc_ = false;
     bool isScaCoeffsCalc_ = false;