Merge branch 'master' of github.com:ovidiopr/scattnlay

README.md

@@ -59,6 +59,7 @@ Compilation options
  - **make deb** - Generate a deb package for Python extension
  - **make standalone** - Create standalone programs (scattnlay and fieldnlay)
  - **make clean** - Delete temporal files
+  
 
 Binary install:
 --------------
@@ -79,6 +80,8 @@ You can also install it from PyPi via
 sudo pip install python-scattnlay
 ```
 
+You can also ```git clone``` and ```pip install -e .``` to develop python package.
+
 Use:
 ----
 

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

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of python-scattnlay
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite at least one of the following references:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+#        calculation of electromagnetic near-field for a multilayered
+#        sphere," Computer Physics Communications, vol. 214, May 2017,
+#        pp. 225-230.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http:#www.gnu.org/licenses/>.
+
+import scattnlay
+from scattnlay import fieldnlay
+from scattnlay import scattnlay
+import numpy as np
+from matplotlib import pyplot as plt
+
+npts = 151
+factor = 3. # plot extent compared to sphere radius
+index_H2O = 1.33+0.j
+
+WL = 0.532 #mkm
+total_r = 3 #mkm
+
+nm = 1.0 # host medium
+x = 2.0 * np.pi * np.array([total_r / 4.0 * 3.0, total_r], dtype=np.float64) / WL
+m = np.array((index_H2O, index_H2O), dtype=np.complex128) / nm
+
+print("x =", x)
+print("m =", m)
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+    np.array([x]), np.array([m]))
+print("Qsca = " + str(Qsca))
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+    np.array([x]), np.array([m]), mp=True)
+print("mp Qsca = " + str(Qsca))
+
+scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
+zero = np.zeros(npts*npts, dtype = np.float64)
+
+coordX, coordZ = np.meshgrid(scan, scan)
+coordX.resize(npts * npts)
+coordZ.resize(npts * npts)
+coordY = zero
+
+terms, E, H = fieldnlay(
+    np.array([x]), np.array([m]),
+    coordX, coordY, coordZ,
+    mp=True
+)
+Ec = E[0, :, :]
+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$'
+plt.imshow(Eabs_data,
+           cmap='jet',
+           vmin=0., vmax=14
+
+           )
+print(np.min(Eabs_data), np.max(Eabs_data))
+plt.savefig("R"+str(total_r)+"mkm_mp.jpg")
+# plt.show()

examples/field-Ag-flow.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
@@ -26,7 +26,7 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-# This test case calculates the electric field in the 
+# This test case calculates the electric field in the
 # E-k plane, for an spherical Ag nanoparticle.
 
 from scattnlay import fieldnlay, scattnlay
@@ -34,29 +34,30 @@ from fieldplot import fieldplot
 
 import numpy as np
 import cmath
+
 # # a)
-#WL=400 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.0 + 10.0j
+# WL=400 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.0 + 10.0j
 
 # # b)
-#WL=400 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.0 + 1.0j
-
-# # c)
-# WL=354 #nm
+# WL=400 #nm
 # core_r = WL/20.0
-# epsilon_Ag = -2.0 + 0.28j
+# epsilon_Ag = -2.0 + 1.0j
+
+# c)
+WL = 354  # nm
+core_r = WL / 20.0
+epsilon_Ag = -2.0 + 0.28j
 
 # d)
-#WL=367 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.71 + 0.25j
+# WL=367 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.71 + 0.25j
 
-WL=500 #nm
-core_r = 50.0
-epsilon_Ag = 4.0 
+# WL=500 #nm
+# core_r = 615.0
+# epsilon_Ag = 4.0
 
 
 index_Ag = np.sqrt(epsilon_Ag)
@@ -64,36 +65,37 @@ index_Ag = np.sqrt(epsilon_Ag)
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
-x = 2.0*np.pi*np.array([core_r/4.0*3.0, core_r], dtype = np.float64)/WL
+x = 2.0 * np.pi * np.array([core_r / 4.0 * 3.0, core_r], dtype=np.float64) / WL
 
-m = np.array((index_Ag, index_Ag), dtype = np.complex128)/nm
+m = np.array((index_Ag, index_Ag), dtype=np.complex128) / nm
 
-print "x =", x
-print "m =", m
+print("x =", x)
+print("m =", m)
 
-comment='bulk-Ag-flow'
-WL_units='nm'
-npts = 151
-factor=2.1
+comment = 'bulk-WL' + str(WL) + 'nm_r' + str(core_r) + 'nm_epsilon' + str(epsilon_Ag) + '-flow'
+WL_units = 'nm'
+npts = 251
+factor = 2.1
 flow_total = 9
-#flow_total = 21
-#flow_total = 0
-#crossplane='XZ'
-#crossplane='YZ'
-crossplane='XY'
+# flow_total = 21
+# flow_total = 0
+crossplane = 'XZ'
+# crossplane='YZ'
+# crossplane='XY'
 
 # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-field_to_plot='Eabs'
-#field_to_plot='angleEx'
+field_to_plot = 'Pabs'
+# field_to_plot='angleEx'
 
 
 import matplotlib.pyplot as plt
-fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+
+fig, axs = plt.subplots(1, 1)  # , sharey=True, sharex=True)
 fig.tight_layout()
-fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-          subplot_label=' ',is_flow_extend=False)
+fieldplot(fig, axs, x, m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+          subplot_label=' ', is_flow_extend=False)
 
-#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+# fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
 
 # for ax in axs:
 #     ax.locator_params(axis='x',nbins=5)
@@ -101,8 +103,8 @@ fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts,
 
 fig.subplots_adjust(hspace=0.3, wspace=-0.1)
 
-plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
-                    +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+plt.savefig(comment + "-R" + str(int(round(x[-1] * WL / 2.0 / np.pi))) + "-" + crossplane + "-"
+            + field_to_plot + ".pdf", pad_inches=0.02, bbox_inches='tight')
 
 plt.draw()
 

examples/field-SiAgSi-flow.py

@@ -26,7 +26,7 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-# This test case calculates the electric field in the 
+# This test case calculates the electric field in the
 # E-k plane, for an spherical Si-Ag-Si nanoparticle. Core radius is 17.74 nm,
 # inner layer 23.31nm, outer layer 22.95nm. Working wavelength is 800nm, we use
 # silicon epsilon=13.64+i0.047, silver epsilon= -28.05+i1.525
@@ -72,11 +72,19 @@ outer_r = inner_r+outer_width
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
+# WL=354 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.0 + 0.28j
+# index_Ag = np.sqrt(epsilon_Ag)
+# x = 2.0*np.pi*np.array([core_r/3., core_r/2., core_r], dtype = np.float64)/WL
+# m = np.array((index_Ag, index_Ag, index_Ag), dtype = np.complex128)/nm
+
+
 x = 2.0*np.pi*np.array([core_r, inner_r, outer_r], dtype = np.float64)/WL
 m = np.array((index_Si, index_Ag, index_Si), dtype = np.complex128)/nm
 
-print "x =", x
-print "m =", m
+print("x =", x)
+print("m =", m)
 
 npts = 501
 factor=2.2

examples/fieldplot.py

@@ -103,7 +103,7 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
             Ec, Hc = E[0, 0, :], H[0, 0, :]
             Eth = max(np.absolute(Ec)) / 1e10
             Hth = max(np.absolute(Hc)) / 1e10
-            for i in xrange(0, len(Ec)):
+            for i in range(0, len(Ec)):
                 if abs(Ec[i]) < Eth:
                     Ec[i] = 0 + 0j
                 if abs(Hc[i]) < Hth:
@@ -186,10 +186,12 @@ def GetField(crossplane, npts, factor, x, m, pl):
     terms, E, H = fieldnlay(np.array([x]), np.array([m]), coordX, coordY, coordZ, pl=pl)
     Ec = E[0, :, :]
     Hc = H[0, :, :]
-    P = []
-    P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real,
-                     range(0, len(E[0]))))
-
+    P = np.array(list(map(lambda n: np.linalg.norm(np.cross(Ec[n],
+                                                            np.conjugate(Hc[n])
+                                                            # Hc[n]
+                                                            )).real,
+                     range(0, len(E[0])))))
+    print(P)
     # for n in range(0, len(E[0])):
     #     P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
     return Ec, Hc, P, coordPlot1, coordPlot2
@@ -199,7 +201,9 @@ def GetField(crossplane, npts, factor, x, m, pl):
 def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
               field_to_plot='Pabs', npts=101, factor=2.1, flow_total=11,
               is_flow_extend=True, pl=-1, outline_width=1, subplot_label=' '):
-
+    # print(fig, ax, x, m, WL, comment, WL_units, crossplane,
+    #       field_to_plot, npts, factor, flow_total,
+    #       is_flow_extend, pl, outline_width, subplot_label)
     Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl)
     Er = np.absolute(Ec)
     Hr = np.absolute(Hc)
@@ -209,21 +213,22 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
 
         if field_to_plot == 'Pabs':
             Eabs_data = np.resize(P, (npts, npts)).T
-            label = r'$\operatorname{Re}(E \times H)$'
+            label = r'$\operatorname{Re}(E \times H^*)$'
         elif field_to_plot == 'Eabs':
-            # Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
-            # label = r'$|E|$'
+            Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
+            label = r'$|E|$'
             # Eabs = np.real(Hc[:, 0])
             # label = r'$Re(H_x)$'
             # Eabs = np.real(Hc[:, 1])
             # label = r'$Re(H_y)$'
-            Eabs = np.real(Ec[:, 1])
-            label = r'$Re(E_y)$'
+            # Eabs = np.real(Ec[:, 1])
+            # label = r'$Re(E_y)$'
             # Eabs = np.real(Ec[:, 0])
             # label = r'$Re(E_x)$'
-            Eabs_data = np.resize(Eabs, (npts, npts))
+            Eabs_data = np.resize(Eabs, (npts, npts)).T
         elif field_to_plot == 'Habs':
             Habs = np.sqrt(Hr[:, 0]**2 + Hr[:, 1]**2 + Hr[:, 2]**2)
+            Habs = 376.730313667 * Habs # scale to free space impedance
             Eabs_data = np.resize(Habs, (npts, npts)).T
             label = r'$|H|$'
         elif field_to_plot == 'angleEx':

setup.py

@@ -64,12 +64,12 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
       ext_modules = [Extension("scattnlay_dp",
                                ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
                                language = "c++",
-                               include_dirs = [np.get_include(), pb.get_include()], 
+                               include_dirs = [np.get_include(), pb.get_include()],
                                extra_compile_args=['-std=c++11']),
                      Extension("scattnlay_mp",
                                ["src/nmie.cc", "src/nmie-pybind11.cc", "src/pb11_wrapper.cc"],
                                language = "c++",
-                               include_dirs = [np.get_include(), pb.get_include()], 
+                               include_dirs = [np.get_include(), pb.get_include()],
                                extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])]
 )
 

src/nmie-impl.cc

@@ -58,7 +58,7 @@
 #include <stdexcept>
 #include <vector>
 
-namespace nmie {  
+namespace nmie {
   //helper functions
 
 
@@ -73,11 +73,11 @@ namespace nmie {
     //return x >= 0 ? (x + 0.5).convert_to<int>():(x - 0.5).convert_to<int>();
   }
   template<typename T>
-  inline std::complex<T> my_exp(const std::complex<T>& x) { 
+  inline std::complex<T> my_exp(const std::complex<T>& x) {
     using std::exp; // use ADL
     T const& r = exp(x.real());
-    return std::polar(r, x.imag()); 
-  }  
+    return std::polar(r, x.imag());
+  }
 
 
   //class implementation
@@ -455,7 +455,7 @@ namespace nmie {
   void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
                                   std::vector<std::complex<FloatType> >& Psi,
                                   std::vector<std::complex<FloatType> >& Zeta) {
-  
+
     std::complex<FloatType> c_i(0.0, 1.0);
     std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1);
 
@@ -528,7 +528,7 @@ namespace nmie {
   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> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n,
                                     std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n) {
 
     // using eq 4.50 in BH
@@ -995,7 +995,7 @@ namespace nmie {
       E[i] = c_zero;
       H[i] = c_zero;
     }
-    
+
     if (Rho > size_param_.back()) {
       l = size_param_.size();
       ml = c_one;
@@ -1098,12 +1098,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;