Browse Source

correct dielectiric sphere field plot

Konstantin Ladutenko 10 years ago
parent
commit
f786a93665
3 changed files with 37 additions and 21 deletions
  1. 24 9
      compare.cc
  2. 7 7
      go.sh
  3. 6 5
      tests/python/field-dielectric-sphere.py

+ 24 - 9
compare.cc

@@ -234,22 +234,26 @@ int main(int argc, char *argv[]) {
       }
     }
     // Field testing
-    double from_coord = -3.0, to_coord = 3000.0;
     //double size=2.0*PI*1.0/6.0;
     double size=0.001;
     double R = size/(2.0*PI);
+    double from_coord = -3.0*size, to_coord = 3.0*size;
     std::vector<double> range;
-    // for (int i = 0; i < samples; ++i) {
-    //range.push_back( from_coord + (to_coord-from_coord)/(static_cast<double>(samples)-1)*i );
+    int samples = 1251;
+    for (int i = 0; i < samples; ++i) {
+    range.push_back( from_coord + (to_coord-from_coord)/(static_cast<double>(samples)-1)*i );
     //range.push_back(size*0.01);
-    range.push_back(size*0.99999);
+    //range.push_back(size*0.99999);
     //range.push_back(R/2.0);
-    range.push_back(size*1.00001);
+    //range.push_back(size*1.00001);
     //range.push_back(3);
     //printf("r=%g  ", range.back());
-    //}
+    }
+    // range.push_back(size*0.99999999);
+    // range.push_back(R/2.0);
+    // range.push_back(size*1.00000001);
     //printf("r/2 = %g\n", R/2.0);
-    int samples = range.size();
+    //int samples = range.size();
     std::vector<double> zero(samples, 0.0);
     std::vector<double> Xp, Yp, Zp;
      // X line
@@ -282,14 +286,25 @@ int main(int argc, char *argv[]) {
     nmie::nField( L,  pl,  x,  m, nmax,  ncoord,  Xp,  Yp,  Zp, E, H);
     double sum_e = 0.0, sum_h = 0.0;
     printf ("Field total sum ()\n");
+    double min_E, max_E;
+    for (auto c:E[0]) {
+      sum_e+=std::abs(pow2(c));
+    }
+    min_E = sum_e;
+    max_E = sum_e;
+
     for (auto f:E) {
       sum_e = 0.0;
       for (auto c:f) {
 	sum_e+=std::abs(pow2(c));
-	printf("component: %g + %g i\n", std::real(c), std::imag(c));
+	//printf("component: %g + %g i\n", std::real(c), std::imag(c));
       }
-      printf("Field E=%g\n", std::sqrt(std::abs(sum_e)));
+      if (sum_e > max_E) max_E = sum_e;
+      if (sum_e < min_E) min_E = sum_e;
+    
+      //printf("Field E=%g\n", std::sqrt(std::abs(sum_e)));
     }
+    printf("Min E = %g; max E =%g", min_E, max_E);
     // for (auto f:H) {
     //   sum_h = 0.0;
     //   for (auto c:f) sum_h+=std::abs(pow2(c));

+ 7 - 7
go.sh

@@ -94,13 +94,13 @@ time  ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4  $PROGRAM -l 5 0.4642 1.
 # ./$file.bin
 
 ### Cython
-rm scattnlay.so
-export CFLAGS='-std=c++11'
-python setup.py build_ext --inplace
-cp scattnlay.so tests/python/
-cd tests/python/
-./lfield.py
-# ./field-dielectric-sphere..py
+# rm scattnlay.so
+# export CFLAGS='-std=c++11'
+# python setup.py build_ext --inplace
+# cp scattnlay.so tests/python/
+# cd tests/python/
+# ./lfield.py
+# ./field-dielectric-sphere.py
 # ./field.py
 # ./test01.py
 # ./test01-wrapper.py

+ 6 - 5
tests/python/field-dielectric-sphere.py

@@ -47,7 +47,7 @@ x[0, 0] = 0.001
 #x[0, 1] = 2.0*np.pi*0.06/1.064
 
 m = np.ones((1, 1), dtype = np.complex128)
-m[0, 0] = 4.0
+m[0, 0] = 2.0
 #m[0, 1] = (0.565838 + 7.23262j)/1.3205
 
 npts = 501
@@ -63,12 +63,13 @@ coord = np.vstack((coordX, coordY, coordZ)).transpose()
 
 terms, E, H = fieldnlay(x, m, coord)
 
-Er = np.absolute(E)
+Er = E
 
 # |E|/|Eo|
-Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
-Eh = Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2
+#Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+Eh = np.absolute(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
 
+print("max: "+str(Eh)+" min: "+str(min(E)))
 result = np.vstack((coordX, coordY, coordZ, Eh)).transpose()
 
 try:
@@ -106,7 +107,7 @@ try:
     cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
     cbar.ax.set_yticklabels(['%3.1e' % (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)
+    fig.text(pos[0] - 0.02, 0.925, 'E$^2$/E$_0$$^2$', fontsize = 14)
 
     plt.xlabel('X')
     plt.ylabel('Y')