Coulomb works fine with Maxwell stress tensor

examples/calc-expansion-coeffs-spectra.py

@@ -83,6 +83,7 @@ plt.xlabel("D, nm")
 plt.ylabel("Mie coefficients")
 plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+ext,
             pad_inches=0.02, bbox_inches='tight')
+plt.show()
 plt.clf()
 plt.close()
 

examples/example-eval-force.cc

@@ -85,8 +85,8 @@ int main(int argc, char *argv[]) {
         // std::cout<<"\tforce:\t" <<F[0]<<"\t"<< F[1] <<"\t"<<F[2] << std::endl;
 
         // std::cout<<"clean"<<std::endl;
-        auto F = shell.IntegrateGauss(2.54,12.03);
-        //std::cout<<"\tcharge:\t" <<F<< std::endl;
+        auto F = shell.IntegrateGaussSimple(2.54,2.03);
+        std::cout<<"\tcharge:\t" <<F<< std::endl;
         
         // auto F1 = shell.IntegrateByComp();
         // std::cout<<"F: " <<F1[0]<<", "<< F1[1] <<", "<<F1[2] << std::endl;        

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

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2018  Konstantin Ladutenko <kostyfisik@gmail.com>
@@ -50,31 +50,68 @@ WL = 550
 x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
 m[0] = index_Ag/nm
 
-R = x[0]*1.31
+R_st = x[0]*4.11
+#R_st = 0.31
+
+dx = R_st*4.0
+charge = 1.0
 
 comment='bulk-NP-WL'+str(WL)+WL_units
 
-quad_ord = 19
-#quad_ord = 131
+#quad_ord = 3
+# quad_ord = 19
+quad_ord = 131
 
-coord = quadpy.sphere.Lebedev(3).points
+#coord = quadpy.sphere.Lebedev(3).points
 # coord = np.vstack((coordX, coordY, coordZ)).transpose()
+def field(coord):
+    E = []
+    for rr in coord:
+        shift = np.array([dx, 0.0, 0.0])
+        unit = (rr-shift)/np.linalg.norm(rr-shift)
+        norm = np.linalg.norm(rr-shift)
+        amp = charge/(4*np.pi*(norm**2))
+        Eloc = amp*unit
+        
+        shift = np.array([0.0, 0.0, 0.0])
+        unit = (rr-shift)/np.linalg.norm(rr-shift)
+        norm = np.linalg.norm(rr-shift)
+        amp = charge/(4*np.pi*(norm**2))
+        Eloc += amp*unit
+        
+        E.append(Eloc)
+    E = np.array(E)
+    return E
+
+def gauss(in_coord):
+    coord = in_coord.T
+    E_all = field(coord)
+    unit_all = coord/ R
+    P = np.array([
+        np.dot(E,unit)
+        for unit,E in zip(unit_all,E_all)
+        ])
+    return P.T
+
+def dipole(coord):
+    H = np.array([0.0,0.0,0.0]*len(coord))
+    E = field(coord)
+    return E,H
+
 def force(in_coord):
     coord = in_coord.T
-    terms, Eall, Hall = fieldnlay(np.array([x]), np.array([m]), coord)
+    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, :, :]
-      # std::vector<double> P = (1/(2.0))
-      #   *real(
-      #         dot(unit,E)*vconj(E) +
-      #         dot(unit,H)*vconj(H) +
-      #         (-1.0/2.0)*(dot(E,vconj(E))
-      #                     +dot(H,vconj(H))
-      #                     )*unit
-      #         );
+
+    E_all, H_all = dipole(coord)
+    # print(coord)
+    # print(E_all)
+    # print(H_all)
+
     unit_all = coord/ R
     P = np.array([
-        ( (1/(2.0))
+        ( (1.0/(2.0)*(4.0*np.pi))
         *np.real(
             np.dot(unit,E)*np.conj(E) +
             np.dot(unit,H)*np.conj(H) +
@@ -94,7 +131,7 @@ def poynting(in_coord):
     H_all = Hall[0, :, :]
     unit_all = coord/ R
     P = np.array([
-        ( ( 1/(2.0) )
+        ( ( 1.0/(2.0) )
             *np.real(
                 np.cross(E,np.conj(H))
             )
@@ -103,23 +140,40 @@ def poynting(in_coord):
         ])
     return P.T
 
+
 # P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real,
 #                      range(0, len(E[0]))))
-
+R=R_st
 val = quadpy.sphere.integrate(
-#    force
-    poynting
+    force
+#    poynting
     ,
     [0.0, 0.0, 0.0], R,
     quadpy.sphere.Lebedev(quad_ord)
     )
 print(val)
 print("Random increase of integraion sphere radius...")
+R=R_st*3.0
 val = quadpy.sphere.integrate(
-    #force
-    poynting
+    force
+   # poynting
     ,
-    [0.0, 0.0, 0.0], R*2.718,
+    [0.0, 0.0, 0.0], R,
     quadpy.sphere.Lebedev(quad_ord)
     )
 print(val)
+
+R=R_st
+print("\n\nCheck Gauss law")
+val = quadpy.sphere.integrate(gauss,
+    [0.0, 0.0, 0.0], R,
+    quadpy.sphere.Lebedev(quad_ord)
+    )
+print("Charge: ",val)
+print("Random increase of integraion sphere radius...")
+R=R_st*3
+val = quadpy.sphere.integrate(gauss,
+        [0.0, 0.0, 0.0], R,
+    quadpy.sphere.Lebedev(quad_ord)
+    )
+print("Charge: ",val)

src/shell-generator.cc

@@ -246,7 +246,7 @@ namespace shell_generator {
       auto vert = face_centers_[i];
       auto E0 = field(charge, shift, vert);
       // std::cout << "E0: ";
-      for (auto component : E0) std::cout << component << " ";
+      // for (auto component : E0) std::cout << component << " ";
       // std::cout << std::endl;
       // Vector to unit product
       double r = norm(vert);