problem field

examples/example-eval-force.cc

@@ -66,7 +66,7 @@ int main(int argc, char *argv[]) {
     //shell.Refine();
 
     //double scale = 2.0*pi*(core_width + inner_width +  outer_width)/WL*8.0001;  //Integration sphere radius.
-    double scale = 2.0*pi*(outer_width)/WL*1.0001;  //Integration sphere radius.
+    double scale = 2.0*pi*(outer_width)/WL*1.4001;  //Integration sphere radius.
     //double scale = 1.0001;  //Integration sphere radius.
     
     shell.Rescale(scale);
@@ -80,10 +80,10 @@ int main(int argc, char *argv[]) {
     auto E = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldE());
     auto H = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldH());
     shell.SetField(E,H);
-    auto F = shell.Integrate();
-    std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
-    F = shell.IntegrateByComp();
-    std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl;
+    // auto F = shell.Integrate();
+    // std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
+    auto F = shell.IntegrateByComp();
+    //std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl;
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;

examples/field-Ag-flow.py

@@ -46,16 +46,21 @@ 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 
+
+
 index_Ag = np.sqrt(epsilon_Ag)
 # n1 = 1.53413
 # n2 = 0.565838 + 7.23262j
@@ -74,17 +79,17 @@ print "m =", m
 
 comment='bulk-Ag-flow'
 WL_units='nm'
-npts = 501
+npts = 151
 factor=2.1
-flow_total = 39
+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='Pabs'
+field_to_plot='Eabs'
 #field_to_plot='angleEx'
 
 

examples/fieldplot.py

@@ -212,9 +212,17 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
             Eabs_data = np.resize(P, (npts, npts)).T
             label = r'$\operatorname{Re}(E \times H)$'
         elif field_to_plot == 'Eabs':
-            Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
-            Eabs_data = np.resize(Eabs, (npts, npts)).T
-            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[:, 0])
+            # label = r'$Re(E_x)$'
+            Eabs_data = np.resize(Eabs, (npts, npts))
         elif field_to_plot == 'Habs':
             Habs = np.sqrt(Hr[:, 0]**2 + Hr[:, 1]**2 + Hr[:, 2]**2)
             Eabs_data = np.resize(Habs, (npts, npts)).T
@@ -263,7 +271,7 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         if 'angle' in field_to_plot:
             cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
         else:
-            cbar.ax.set_yticklabels(['%3.2f' % (a) for a in scale_ticks])
+            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
@@ -278,8 +286,8 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
             ax.set_xlabel(r'$Z,\lambda$'+WL_units)
             ax.set_ylabel(r'$Y:X,\lambda$'+WL_units)
         elif crossplane == 'XY':
-            ax.set_xlabel('Y, ' + WL_units, labelpad=lp1)
-            ax.set_ylabel('X, ' + WL_units, labelpad=lp2)
+            ax.set_xlabel('X, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('Y, ' + WL_units, labelpad=lp2)
         # # This part draws the nanoshell
         from matplotlib import patches
         from matplotlib.path import Path
@@ -367,5 +375,5 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
     finally:
         terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
             np.array([x]), np.array([m]))
-        print("Qabs = " + str(Qabs))
+        print("Qsca = " + str(Qsca))
     #

src/shell-generator.cc

@@ -259,10 +259,12 @@ namespace shell_generator {
       auto vert = vertices_[i];
       std::cout <<"E "<<E[0]<<", "<< E[1] <<", "<<E[2] << std::endl;
       std::cout <<"H "<<H[0]<<", "<< H[1] <<", "<<H[2] << std::endl;
-      std::cout <<"vert "<<vert[0]<<", "<< vert[1] <<", "<<vert[2] << std::endl<<std::endl;
+      std::cout <<"vert "<<vert[0]<<", "<< vert[1] <<", "<<vert[2] << std::endl;
       // Vector to unit product
       double r = norm(vert);
       std::vector<double> n = { vert[0]/r, vert[1]/r, vert[2]/r};
+      //std::cout<<"norm: " <<n[0]<<", "<< n[1] <<", "<<n[2] << std::endl;
+          
       // std::cout << norm(unit) << std::endl;
       //const double pi = 3.1415926535897932384626433832795;
       const std::vector< std::vector<double> > d = {{1.0, 0.0, 0.0},
@@ -272,18 +274,19 @@ namespace shell_generator {
       std::vector<double> F = {0.0, 0.0, 0.0};
       std::complex<double> S(0.0);
       for (int ii = 0; ii < 3; ++ii)
-        S += E[ii]*E[ii] + H[ii]*H[ii];
+        S += E[ii]*std::conj(E[ii]) + H[ii]*std::conj(H[ii]);
       std::vector< std::vector<double> >  T = {{0.0, 0.0, 0.0},
                                {0.0, 0.0, 0.0},
                                {0.0, 0.0, 0.0}};
       for (int i = 0; i < 3; ++i) {
         for (int j = 0; j < 3; ++j) {
           T[i][j] = (1/(2.0/* *4.0*pi */))
-            * real(E[i]*E[j] + H[i]*H[j]
+            * real(E[i]*std::conj(E[j]) + H[i]*std::conj(H[j])
                    -1.0/2.0*S*d[i][j]);
           F[i] += T[i][j]*n[j];
         }
       }
+      std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
       integral = integral + per_vertice_area_*F;
     }
     return integral;
@@ -485,7 +488,7 @@ namespace shell_generator {
 
     //std::vector< std::vector<double> > points_debug = {{1,0,0},{-1,0,0}};
     //std::vector< std::vector<double> > points_debug = {{0,1,0},{0,-1,0}};
-    std::vector< std::vector<double> > points_debug = {{1,1,0},{-1,-1,0},{-1,1,0},{1,-1,0}};
+    std::vector< std::vector<double> > points_debug = {{1,1,0},{1,-1,0},{-1,-1,0},{-1,1,0}};
     //std::vector< std::vector<double> > points_debug = {};
     // std::vector< std::vector<double> > points_debug = {{0,0,1},{0,0,-1}};
     vertices_ = std::move(points_debug);