small update

examples/calc-SiAgSi-Qabs.py

@@ -93,6 +93,24 @@ def SetXM(design):
         # m[0] = index_Ag
         # m[1] = index_Ag
         return x, m, WL
+    elif design==6:
+        WL=1052 #nm
+        core_r = 140.0
+        #core_r = 190.0
+        core_r = 204.2
+        epsilon_Si = 12.7294053067+0.000835315166667j
+        index_Si = np.sqrt(epsilon_Si)
+        x = np.ones((1), dtype = np.float64)
+        x[0] = 2.0*np.pi*core_r/WL
+        m = np.ones((1), dtype = np.complex128)
+        m[0] = index_Si
+        # x = np.ones((2), dtype = np.float64)
+        # x[0] = 2.0*np.pi*core_r/WL/4.0*3.0
+        # x[1] = 2.0*np.pi*core_r/WL
+        # m = np.ones((2), dtype = np.complex128)
+        # m[0] = index_Ag
+        # m[1] = index_Ag
+        return x, m, WL
 
 
     core_r = core_width
@@ -125,29 +143,49 @@ def SetXM(design):
 #design = 1 #AgSi
 #design = 2
 #design = 3
-design = 4   # WL=800
-comment='SiAgSi-flow'
+# design = 4   # WL=800
+# comment='SiAgSi-flow'
 #design = 5   # Bulk Ag
 # comment='bulk-Ag-flow'
+design = 6   # WL=800
+comment='Si-flow'
 x, m, WL = SetXM(design)
 
 WL_units='nm'
-print "x =", x
+print "x =", x[-1]
 print "m =", m
-npts = 101
+npts = 501
 factor=2.1
 flow_total = 39
 #flow_total = 21
 #flow_total = 0
-crossplane='XZ'
+#crossplane='XZ'
+crossplane='XYZ'
 #crossplane='YZ'
 #crossplane='XY'
 
 # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-field_to_plot='Pabs'
+field_to_plot='Eabs'
 #field_to_plot='angleEx'
+#field_to_plot='Pabs'
+
+import matplotlib.pyplot as plt
+fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+fig.tight_layout()
+fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+          subplot_label=' ',is_flow_extend=False)
+
+fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+
+plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
+                    +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+
+plt.draw()
+
+#    plt.show()
 
-fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total)
+plt.clf()
+plt.close()
 
 
 

examples/fieldplot.py

@@ -200,6 +200,7 @@ 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 (x,m)
     Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl)
     Er = np.absolute(Ec)
     Hr = np.absolute(Hc)

examples/read-spectra.h

@@ -0,0 +1,46 @@
+#ifndef SRC_READ_SPECTRA_READ_SPECTRA_H_
+#define SRC_READ_SPECTRA_READ_SPECTRA_H_
+/**
+ * @file   read-spectra.h
+ * @author Konstantin Ladutenko <kostyfisik at gmail (.) com>
+ * @date   Wed Mar 11 11:19:34 2015
+ * @copyright 2015 Konstantin Ladutenko
+ *
+ * @brief  Read complex spectra from file in format 'WL real imag'
+ * 
+ * read-spectra is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * read-spectra is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with read-spectra.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ */
+#include <vector>
+#include <string>
+#include <complex>
+#include <utility>
+namespace read_spectra {
+  class ReadSpectra {  // will throw for any error
+   public:
+    ReadSpectra& ReadFromFile(std::string filename);
+    ReadSpectra& ResizeToComplex(double from_wl, double to_wl, int samples);
+    ReadSpectra& ToIndex();
+    std::complex<double> at(double wavelength);
+    void PrintData();    
+    std::vector< std::pair< double, std::complex<double> > >&
+      GetIndex(){return data_complex_index_;};
+  private:
+    std::vector< std::vector<double> > data_;
+    std::vector< std::pair< double, std::complex<double> > > data_complex_;
+    std::vector< std::pair< double, std::complex<double> > > data_complex_index_;
+    void PermittivityToIndex();
+  };  // end of class ReadSpectra
+}  // end of namespase read_spectra
+#endif  // SRC_READ_SPECTRA_READ_SPECTRA_H_

src/nmie.cc

@@ -141,11 +141,11 @@ namespace nmie {
 
       ml_mie.RunMieCalculation();
 
-      std::cout
-      	<< std::setprecision(std::numeric_limits<FloatType>::digits10)
-      	<< "Qext = "
-      	<< ml_mie.GetQext()
-      	<< std::endl;
+      // std::cout
+      // 	<< std::setprecision(std::numeric_limits<FloatType>::digits10)
+      // 	<< "Qext = "
+      // 	<< ml_mie.GetQext()
+      // 	<< std::endl;
       
       *Qext = static_cast<double>(ml_mie.GetQext());
       *Qsca = static_cast<double>(ml_mie.GetQsca());

tests/c++/go-speed-test.sh

@@ -10,7 +10,7 @@ rm -f $PROGRAM
 # g++ -Ofast -std=c++11 $file ../../src/nmie.cc  -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 file=speed-test.cc
-#g++ -Ofast -std=c++11 $file ../../src/nmie.cc -DMULTI_PRECISION=20  -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+# g++ -Ofast -std=c++11 $file ../../src/nmie.cc -DMULTI_PRECISION=200  -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 g++ -Ofast -std=c++11 $file ../../src/nmie.cc  -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 #file=speed-test-applied.cc

tests/c++/speed-test.cc

@@ -185,37 +185,38 @@ int main(int argc, char *argv[]) {
     long cpptime_sec, ctime_sec;
     long repeats = 150;
     //HeapProfilerStart("heapprof");    
+    printf("Start\n");
     do {
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
       for (int i = 0; i<repeats; ++i) {
     	nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw,
     			   &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
-	break;
+	// break;
       }
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
       cpptime_nsec = diff(time1,time2).tv_nsec;
       cpptime_sec = diff(time1,time2).tv_sec;
 
-      printf("-- C++ time consumed %lg sec\n", (cpptime_nsec/1e9));
+      printf("-- C++ time consumed %lg sec\n", cpptime_sec+(cpptime_nsec/1e9));
       repeats *= 10;
-      break;
-    } while (cpptime_nsec < 1e8 && ctime_nsec < 1e8);
+      // break;
+    } while (cpptime_sec < 3 && ctime_sec < 3);
 
-        printf("\n");
+        printf("Finish\n");
     
-    if (has_comment) {
-      printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
-    } else {
-      printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
-    }
+    // if (has_comment) {
+    //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    // } else {
+    //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    // }
     
-    if (nt > 0) {
-      printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+    // if (nt > 0) {
+    //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
       
-      for (i = 0; i < nt; i++) {
-        printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
-      }
-    }
+    //   for (i = 0; i < nt; i++) {
+    //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
+    //   }
+    // }
 
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.