Small changes

compare.cc

@@ -42,7 +42,17 @@
 
 timespec diff(timespec start, timespec end);
 const double PI=3.14159265358979323846;
-  template<class T> inline T pow2(const T value) {return value*value;}
+template<class T> inline T pow2(const T value) {return value*value;}
+
+template <class VectorType, int dimensions> inline
+std::vector<VectorType> CrossProduct(std::vector<VectorType>& a, std::vector<VectorType>& b) {
+  if (a.size() != 3 || b.size() != 3) throw std::invalid_argument("Cross product only for 3D vectors!");
+  std::vector<VectorType> r (3);   
+  r[0] = a[1]*b[2]-a[2]*b[1];
+  r[1] = a[2]*b[0]-a[0]*b[2];
+  r[2] = a[0]*b[1]-a[1]*b[0];
+  return r;
+}
 
 //***********************************************************************************//
 // This is the main function of 'scattnlay', here we read the parameters as          //
@@ -235,7 +245,12 @@ int main(int argc, char *argv[]) {
     }
     // Field testing
     //double size=2.0*PI*1.0/6.0;
-    double size=0.001;
+    double WL=354; //nm
+    double core_r = WL/20.0;
+    double r_x = 2.0*PI*core_r/WL;
+
+
+    double size=r_x;
     double R = size/(2.0*PI);
     double from_coord = -3.0*size, to_coord = 3.0*size;
     std::vector<double> range;
@@ -256,14 +271,14 @@ int main(int argc, char *argv[]) {
     //int samples = range.size();
     std::vector<double> zero(samples, 0.0);
     std::vector<double> Xp, Yp, Zp;
-     // X line
-    Xp.insert(Xp.end(), range.begin(), range.end());
-    Yp.insert(Yp.end(), zero.begin(), zero.end());
-    Zp.insert(Zp.end(), zero.begin(), zero.end());
-    //Y line
-    Xp.insert(Xp.end(), zero.begin(), zero.end());
-    Yp.insert(Yp.end(), range.begin(), range.end());
-    Zp.insert(Zp.end(), zero.begin(), zero.end());
+    //  // X line
+    // Xp.insert(Xp.end(), range.begin(), range.end());
+    // Yp.insert(Yp.end(), zero.begin(), zero.end());
+    // Zp.insert(Zp.end(), zero.begin(), zero.end());
+    // //Y line
+    // Xp.insert(Xp.end(), zero.begin(), zero.end());
+    // Yp.insert(Yp.end(), range.begin(), range.end());
+    // Zp.insert(Zp.end(), zero.begin(), zero.end());
     // Z line
     Xp.insert(Xp.end(), zero.begin(), zero.end());
     Yp.insert(Yp.end(), zero.begin(), zero.end());
@@ -271,7 +286,11 @@ int main(int argc, char *argv[]) {
     int ncoord = Xp.size();
     // Test solid sphere
     x = {size};
-    m = {std::complex<double>(2.000000,0.00)};
+
+    std::complex<double> epsilon_Ag(-2.0, 0.28);
+    m = {std::sqrt(epsilon_Ag)};
+    //m = {std::complex<double>(2.000000,0.00)};
+
     //m = {std::complex<double>(1.414213562, 0.00)};
 
     L = x.size();

go.sh

@@ -99,6 +99,7 @@ export CFLAGS='-std=c++11'
 python setup.py build_ext --inplace
 cp scattnlay.so tests/python/
 cd tests/python/
+#./field-Ag-flow.py
 # ./lfield.py
 # ./field-dielectric-sphere.py
 # ./field.py

nmie.cc

@@ -1406,7 +1406,7 @@ namespace nmie {
       E[i] = c_zero;
       H[i] = c_zero;
     }
-
+    
     if (Rho > size_param_.back()) {
       l = size_param_.size();
       ml = c_one;

tests/python/field-Ag-flow.py

@@ -125,16 +125,13 @@ coordY = np.zeros(npts*npts, dtype = np.float64)
 
 coord = np.vstack((coordX, coordY, coordZ)).transpose()
 #coord = np.vstack((coordY, coordX, coordZ)).transpose()
+#coord = np.vstack((coordY, coordX, coordZ)).transpose()
+
 
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
 terms, E, H = fieldnlay(x, m, coord)
 Er = np.absolute(E)
 Hr = np.absolute(H)
-P=[]
-for n in range(0, len(E[0])):
-    P.append(np.linalg.norm( np.cross(E[0][n], np.conjugate(H[0][n]) ).real/2 ))
-
-print(min(P))
 
 # |E|/|Eo|
 Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
@@ -142,6 +139,12 @@ Ec = E[0, :, :]
 Hc = H[0, :, :]
 Eangle = np.angle(E[0, :, 0])/np.pi*180
 
+P=[]
+for n in range(0, len(E[0])):
+    P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
+
+print(min(P))
+
 Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
 Hangle = np.angle(H[0, :, 1])/np.pi*180
 
@@ -169,10 +172,10 @@ try:
     scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi/nm, max(coordZ)*WL/2.0/np.pi/nm, npts)
 
     # Define scale ticks
-    # min_tick = min(min_tick, np.amin(Eabs_data))
-    # max_tick = max(max_tick, np.amax(Eabs_data))
+    min_tick = np.amin(Eabs_data)
+    max_tick = np.amax(Eabs_data)
     # scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
-    #scale_ticks = np.linspace(min_tick, max_tick, 11)
+    scale_ticks = np.linspace(min_tick, max_tick, 11)
 
     # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
     #ax.set_title('Eabs')
@@ -185,13 +188,13 @@ try:
     ax.axis("image")
 
     # # Add colorbar
-    # cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
-    # cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
+    cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
+    cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
     # pos = list(cbar.ax.get_position().bounds)
     # fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
 
-    plt.xlabel('Z, nm')
-    plt.ylabel('X, nm')
+    # plt.xlabel('Z, nm')
+    # plt.ylabel('X, nm')
 
     # This part draws the nanoshell
     from matplotlib import patches
@@ -201,6 +204,7 @@ try:
 
     from matplotlib.path import Path
     #import matplotlib.patches as patches
+
     flow_total = 31
     for flow in range(0,flow_total):
         flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
@@ -214,6 +218,7 @@ try:
         path = Path(verts, codes)
         patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
         ax.add_patch(patch)
+
     # # Start powerflow lines in the middle of the image
     # flow_total = 131
     # for flow in range(0,flow_total):
@@ -229,7 +234,7 @@ try:
     #     patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
     #     ax.add_patch(patch)
  
-    # plt.savefig("SiAgSi.png")
+    plt.savefig("Ag-flow.png")
     plt.draw()
 
     plt.show()

tests/python/field-SiAgSi-flow.py

@@ -84,6 +84,8 @@ epsilon_Ag = -28.05 + 1.525j
 index_Si = np.sqrt(epsilon_Si)
 index_Ag = np.sqrt(epsilon_Ag)
 
+print(index_Si)
+print(index_Ag)
 # # Values for 800 nm, taken from http://refractiveindex.info/
 # index_Si = 3.69410 + 0.0065435j
 # index_Ag = 0.18599 + 4.9886j
@@ -114,7 +116,7 @@ m[0, 2] = index_Si/nm
 print "x =", x
 print "m =", m
 
-npts = 1281
+npts = 241
 
 factor=2.5
 scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
@@ -124,7 +126,7 @@ coordX.resize(npts*npts)
 coordZ.resize(npts*npts)
 coordY = np.zeros(npts*npts, dtype = np.float64)
 
-coord = np.vstack((coordX, coordX, coordZ)).transpose()
+coord = np.vstack((coordX, coordY, coordZ)).transpose()
 
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
 terms, E, H = fieldnlay(x, m, coord)
@@ -137,6 +139,10 @@ Ec = E[0, :, :]
 Hc = H[0, :, :]
 Eangle = np.angle(E[0, :, 0])/np.pi*180
 
+P=[]
+for n in range(0, len(E[0])):
+    P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
+
 Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
 Hangle = np.angle(H[0, :, 1])/np.pi*180
 
@@ -150,6 +156,7 @@ try:
     min_tick = 0.0
     max_tick = 1.0
 
+    # Eabs_data = np.resize(P, (npts, npts)).T
     Eabs_data = np.resize(Eabs, (npts, npts)).T
     # Eangle_data = np.resize(Eangle, (npts, npts)).T
     # Habs_data = np.resize(Habs, (npts, npts)).T
@@ -164,6 +171,7 @@ try:
     # Define scale ticks
     min_tick = min(min_tick, np.amin(Eabs_data))
     max_tick = max(max_tick, np.amax(Eabs_data))
+    #max_tick = 5
     # scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
     scale_ticks = np.linspace(min_tick, max_tick, 11)
 
@@ -171,7 +179,7 @@ try:
     ax.set_title('Eabs')
     cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = cm.jet,
                     origin = 'lower'
-                    #, vmin = min_tick, vmax = max_tick
+                    , vmin = min_tick, vmax = max_tick
                     , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
                     #,norm = LogNorm()
                     )
@@ -214,7 +222,7 @@ try:
         ax.add_patch(patch)
 
  
-    # plt.savefig("SiAgSi.png")
+    plt.savefig("SiAgSi-flow.png")
     plt.draw()
 
     plt.show()

tests/python/field-SiAgSi.py

@@ -82,7 +82,7 @@ print "m =", m
 
 npts = 281
 
-factor=4.5
+factor=2.5
 scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
 
 coordX, coordZ = np.meshgrid(scan, scan)