small fixes

examples/field-Ag-flow.py

@@ -44,19 +44,19 @@ import cmath
 #core_r = WL/20.0
 #epsilon_Ag = -2.0 + 1.0j
 
-# # c)
-# WL=354 #nm
-# core_r = WL/20.0
-# epsilon_Ag = -2.0 + 0.28j
+# 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=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)
@@ -71,19 +71,19 @@ m = np.array((index_Ag, index_Ag), dtype = np.complex128)/nm
 print( "x =", x)
 print( "m =", m)
 
-comment='bulk-Ag-flow'
+comment='bulk-WL'+str(WL)+'nm_r'+str(core_r)+'nm_epsilon'+str(epsilon_Ag)+'-flow'
 WL_units='nm'
-npts = 151
+npts = 251
 factor=2.1
 flow_total = 9
 #flow_total = 21
 #flow_total = 0
-#crossplane='XZ'
+crossplane='XZ'
 #crossplane='YZ'
-crossplane='XY'
+# crossplane='XY'
 
 # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-field_to_plot='Eabs'
+field_to_plot='Pabs'
 #field_to_plot='angleEx'
 
 

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':

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;