expansion coeffs to python

Makefile

@@ -3,7 +3,7 @@ PYTHON3=`which python3`
 CYTHON=`which cython`
 DESTDIR=/
 PROJECT=python-scattnlay
-VERSION=2.2  # compare with nmie.hpp 
+VERSION=2.2  # compare with nmie.hpp and setup.py
 BUILDIR=$(CURDIR)/debian/$(PROJECT)
 SRCDIR=$(CURDIR)/src
 MULTIPREC=100

examples/example-get-Mie.cc

@@ -25,6 +25,9 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 //   This program returns expansion coefficents of Mie series
+#include "../src/nmie.hpp"
+#include "../src/nmie.hpp"
+#include "../src/nmie-precision.hpp"
 #include <complex>
 #include <cstdio>
 #include <string>
@@ -37,120 +40,126 @@
 int main(int argc, char *argv[]) {
   using namespace nmie ;
   try {
-    read_spectra::ReadSpectra Si_index, Ag_index;
-    read_spectra::ReadSpectra plot_core_index_, plot_TiN_;
-    std::string core_filename("Si-int.txt");
+    //read_spectra::ReadSpectra Si_index;
+    //read_spectra::ReadSpectra plot_core_index_, plot_TiN_;
+   // std::string core_filename("Si-int.txt");
     //std::string core_filename("Ag.txt");
     //std::string TiN_filename("TiN.txt");
-    std::string TiN_filename("Ag-int.txt");
+    //std::string TiN_filename("Ag-int.txt");
     //std::string TiN_filename("Si.txt");
-    std::string shell_filename(core_filename);
+    //std::string shell_filename(core_filename);
 
     nmie::MultiLayerMieApplied<nmie::FloatType> multi_layer_mie;  
-    const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
-    const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
+    const std::complex<double> epsilon_Si(16, 0);
+    //const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
     const std::complex<double> index_Si = std::sqrt(epsilon_Si);
-    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
-    double WL=500; //nm
-    double core_width = 5.27; //nm Si
-    double inner_width = 8.22; //nm Ag
-    double outer_width = 67.91; //nm  Si
-    bool isSiAgSi = true;
-    double delta_width = 25.0;
-    //bool isSiAgSi = false;
-    if (isSiAgSi) {
-      core_width = 5.27; //nm Si
-      inner_width = 8.22; //nm Ag
-      outer_width = 67.91; //nm  Si
-      multi_layer_mie.AddTargetLayer(core_width, index_Si);
-      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
-      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
-    } else {
-      inner_width = 31.93; //nm Ag
-      outer_width = 4.06; //nm  Si
-      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
-      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
-    }
+    //const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
+    double WL=550; //nm
+    //double core_width = 102; //nm Si // radius
+    double core_width = 196/2.0; //nm Si // radius
+    //double inner_width = 8.22; //nm Ag
+    //double outer_width = 67.91; //nm  Si
+    //bool isSiAgSi = true;
+    //double delta_width = 25.0;
+    // //bool isSiAgSi = false;
+    // if (isSiAgSi) {
+      // core_width = 5.27; //nm Si
+      // inner_width = 8.22; //nm Ag
+      // outer_width = 67.91; //nm  Si
+    multi_layer_mie.AddTargetLayer(core_width, index_Si);
+      // multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+      // multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+    // } else {
+    //   inner_width = 31.93; //nm Ag
+    //   outer_width = 4.06; //nm  Si
+    //   multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+    //   multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+    // }
 
-    for (int dd = 0; dd<50; ++dd) {
-      delta_width = dd;
-    FILE *fp;
-    std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat";
-    fp = fopen(fname.c_str(), "w");
+    // for (int dd = 0; dd<50; ++dd) {
+    //   delta_width = dd;
+    // FILE *fp;
+    // std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat";
+    // fp = fopen(fname.c_str(), "w");
 
     multi_layer_mie.SetWavelength(WL);
     multi_layer_mie.RunMieCalculation();
     
-    double Qabs = static_cast<double>(multi_layer_mie.GetQabs());
-    printf("Qabs = %g\n", Qabs);
+    double Qsca = static_cast<double>(multi_layer_mie.GetQsca());
+    printf("Qsca = %g\n", Qsca);//*3.14159*core_width*core_width*1e-6);
     std::vector< std::vector<std::complex<nmie::FloatType> > > aln, bln, cln, dln;
     multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
-    std::vector< std::vector<std::complex<double> > > d_aln =
-      nmie::ConvertComplexVectorVector<double>(aln);
-    std::string str = std::string("#WL ");
-    for (int l = 0; l<d_aln.size(); ++l) {
-      for (int n = 0; n<3; ++n) {
-	str+="|a|^2+|d|^2_ln"+std::to_string(l)+std::to_string(n)+" "
-	  + "|b|^2+|c|^2_ln"+std::to_string(l)+std::to_string(n)+" ";
-      }
-    }
-    str+="\n";
-    fprintf(fp, "%s", str.c_str());
-    fprintf(fp, "# |a|+|d|");
-    str.clear();
-
-    double from_WL = 400;
-    double to_WL = 600;
-    int total_points = 401;
-    double delta_WL = std::abs(to_WL - from_WL)/(total_points-1.0);
-    Si_index.ReadFromFile(core_filename).ResizeToComplex(from_WL, to_WL, total_points)
-      .ToIndex();
-    Ag_index.ReadFromFile(TiN_filename).ResizeToComplex(from_WL, to_WL, total_points)
-      .ToIndex();
-    auto Si_data = Si_index.GetIndex();
-    auto Ag_data = Ag_index.GetIndex();
-    for (int i=0; i < Si_data.size(); ++i) {
-      const double& WL = Si_data[i].first;
-      const std::complex<double>& Si = Si_data[i].second;
-      const std::complex<double>& Ag = Ag_data[i].second;
-      str+=std::to_string(WL);
-      multi_layer_mie.ClearTarget();      
-      if (isSiAgSi) {
-	multi_layer_mie.AddTargetLayer(core_width, Si);
-	multi_layer_mie.AddTargetLayer(inner_width, Ag);
-	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
-      } else {
-	inner_width = 31.93; //nm Ag
-	outer_width = 4.06; //nm  Si
-	multi_layer_mie.AddTargetLayer(inner_width, Ag);
-	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
-      }
-      multi_layer_mie.SetWavelength(WL);
-      multi_layer_mie.RunMieCalculation();
-      multi_layer_mie.GetQabs();
-      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
-      for (int l = 0; l<aln.size(); ++l) {
-	for (int n = 0; n<3; ++n) {
-	  str+=" "+std::to_string(static_cast<double>(pow2(std::abs(aln[l][n]))+
-                                                      pow2(std::abs(dln[l][n]))))
-	    + " "
-	    + std::to_string(static_cast<double>(pow2(std::abs(bln[l][n]))
-                                                 + pow2(std::abs(cln[l][n])) ));
-
-	  // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n]))
-	  // 			  +dln[l][n].real() - pow2(std::abs(dln[l][n])))
-	  //   + " "
-	  //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
-	  // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n])) );
-	}
-      }
-      str+="\n";
-      fprintf(fp, "%s", str.c_str());
-      str.clear();
-    }
 
-    fclose(fp);
-    }
+
+
+
+
+
+ //      nmie::ConvertComplexVectorVector<double>(aln);
+ //    std::string str = std::string("#WL ");
+ //    for (int l = 0; l<d_aln.size(); ++l) {
+ //      for (int n = 0; n<3; ++n) {
+	// str+="|a|^2+|d|^2_ln"+std::to_string(l)+std::to_string(n)+" "
+	//   + "|b|^2+|c|^2_ln"+std::to_string(l)+std::to_string(n)+" ";
+ //      }
+ //    // }
+ //    str+="\n";
+ //    fprintf(fp, "%s", str.c_str());
+ //    fprintf(fp, "# |a|+|d|");
+ //    str.clear();
+
+ //    double from_WL = 400;
+ //    double to_WL = 600;
+ //    int total_points = 401;
+ //    double delta_WL = std::abs(to_WL - from_WL)/(total_points-1.0);
+ //    Si_index.ReadFromFile(core_filename).ResizeToComplex(from_WL, to_WL, total_points)
+ //      .ToIndex();
+ //    Ag_index.ReadFromFile(TiN_filename).ResizeToComplex(from_WL, to_WL, total_points)
+ //      .ToIndex();
+ //    auto Si_data = Si_index.GetIndex();
+ //    auto Ag_data = Ag_index.GetIndex();
+ //    for (int i=0; i < Si_data.size(); ++i) {
+ //      const double& WL = Si_data[i].first;
+ //      const std::complex<double>& Si = Si_data[i].second;
+ //      const std::complex<double>& Ag = Ag_data[i].second;
+ //      str+=std::to_string(WL);
+ //      multi_layer_mie.ClearTarget();      
+ //      if (isSiAgSi) {
+	// multi_layer_mie.AddTargetLayer(core_width, Si);
+	// multi_layer_mie.AddTargetLayer(inner_width, Ag);
+	// multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
+ //      } else {
+	// inner_width = 31.93; //nm Ag
+	// outer_width = 4.06; //nm  Si
+	// multi_layer_mie.AddTargetLayer(inner_width, Ag);
+	// multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
+ //      }
+ //      multi_layer_mie.SetWavelength(WL);
+ //      multi_layer_mie.RunMieCalculation();
+ //      multi_layer_mie.GetQabs();
+ //      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+ //      for (int l = 0; l<aln.size(); ++l) {
+	// for (int n = 0; n<3; ++n) {
+	//   str+=" "+std::to_string(static_cast<double>(pow2(std::abs(aln[l][n]))+
+ //                                                      pow2(std::abs(dln[l][n]))))
+	//     + " "
+	//     + std::to_string(static_cast<double>(pow2(std::abs(bln[l][n]))
+ //                                                 + pow2(std::abs(cln[l][n])) ));
+
+	//   // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n]))
+	//   // 			  +dln[l][n].real() - pow2(std::abs(dln[l][n])))
+	//   //   + " "
+	//   //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
+	//   // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n])) );
+	// }
+ //      }
+ //      str+="\n";
+ //      fprintf(fp, "%s", str.c_str());
+ //      str.clear();
+ //    }
+
+ //    fclose(fp);
+    // }
 
     // WL = 500;
     // multi_layer_mie.SetWavelength(WL);
@@ -158,24 +167,34 @@ int main(int argc, char *argv[]) {
     // multi_layer_mie.GetQabs();
     // multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
 
-    // printf("\n Scattering");
-    // for (int l = 0; l<aln.size(); ++l) {
-    //   int n = 0;
-    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
-    //   n = 1;
-    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
-    //   // n = 2;
-    //   // printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
-    //   // printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
-    //   // printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
-    //   // printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
-    // }
+    printf("\n Scattering");
+    for (int l = 0; l<aln.size(); ++l) {
+      int n = 0;
+      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+      n = 1;
+      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+       n = 2;
+      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+             n = 3;
+      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+             n = 4;
+      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    }
     // printf("\n Absorbtion\n");
     // for (int l = 0; l<aln.size(); ++l) {
     //   if (l == aln.size()-1) printf(" Total ");

examples/fieldplot.py

@@ -132,7 +132,7 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
 
 
 ###############################################################################
-def GetField(crossplane, npts, factor, x, m, pl, mode_n=-1, mode_type=-1):
+def GetField(crossplane, npts, factor, x, m, pl, mode_n=-1, mode_type=-1, inner_only = False):
     """
     crossplane: XZ, YZ, XY, or XYZ (half is XZ, half is YZ)
     npts: number of point in each direction
@@ -189,6 +189,10 @@ def GetField(crossplane, npts, factor, x, m, pl, mode_n=-1, mode_type=-1):
                             )
     Ec = E[0, :, :]
     Hc = H[0, :, :]
+    if inner_only:
+        mask = [ [np.sum(c**2)>x[-1]**2]*3 for c in coord]
+        np.place(Ec, mask, 0)
+        np.place(Hc, mask, 0)
     P = []
     P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real,
                      range(0, len(E[0]))))
@@ -202,9 +206,9 @@ def GetField(crossplane, npts, factor, x, m, pl, mode_n=-1, mode_type=-1):
 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=' ',
-              mode_n=-1, mode_type=-1, isStream = False):
+              mode_n=-1, mode_type=-1, isStream = False,minlength=0.5 , inner_only = False):
     print ("WL=",WL," xm:", x,m)
-    Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl, mode_n, mode_type)
+    Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl, mode_n, mode_type, inner_only)
     Er = np.absolute(Ec)
     Hr = np.absolute(Hc)
     try:
@@ -255,7 +259,7 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
         #scale_ticks = [0.1,0.3,1,3,10, max_tick]
         # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
-        ax.set_title(label)
+        #ax.set_title(label)
         # build a rectangle in axes coords
         ax.annotate(subplot_label, xy=(0.0, 1.1), xycoords='axes fraction',  # fontsize=10,
                     horizontalalignment='left', verticalalignment='top')
@@ -265,22 +269,23 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         #         transform=ax.transAxes)
         cax = ax.imshow(Eabs_data
                         #, interpolation='nearest'
-                        , interpolation='quadric'
+                        , interpolation='none'
+                        #, interpolation='quadric'
                         , cmap=cm.jet,
                         origin='lower', vmin=min_tick, vmax=max_tick, extent=(min(scale_x), max(scale_x), min(scale_z), max(scale_z))
                         # ,norm = LogNorm()
                         )
         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])
+        # # 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
@@ -303,7 +308,7 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         for xx in x:
             r = xx * WL / 2.0 / np.pi
             s1 = patches.Arc((0, 0), 2.0 * r, 2.0 * r,  angle=0.0, zorder=1.8,
-                             theta1=0.0, theta2=360.0, linewidth=outline_width, color='black')
+                             theta1=0.0, theta2=360.0, linewidth=outline_width*1.8, color='black')
             ax.add_patch(s1)
         #
         # for flow in range(0,flow_total):
@@ -394,15 +399,17 @@ def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
         #print(U[:,0])
         X = coordX.reshape((npts, npts))[0,:]*WL/2.0/np.pi
         Z = coordZ.reshape((npts, npts))[:,0]*WL/2.0/np.pi
-        #lw = np.sqrt((V**2+U**2)/(np.max(V)**2+np.max(U)**2))
-        ax.streamplot(X, Z, V, U, color="white", linewidth=0.5,
+        lw = np.sqrt((V**2+U**2)/(np.max(V**2+U**2)))
+        pts2 = npts*3.0
+        
+        ax.streamplot(X, Z, V, U, color="white", linewidth=lw*4.0+1.2,
                       #cmap=plt.cm.inferno,
-                          density=2, arrowstyle='->', arrowsize=1.0)
-
-        #ax.quiver(X[::5], Z[::5], V[::5,::5], U[::5,::5], units='width')
-        print(np.std(field[:,0]),
-              np.std(field[:,1]),
-              np.std(field[:,2])
+                          density=1.1, arrowstyle='->', arrowsize=1.5
+                          ,minlength=minlength
+                          ,maxlength=25
+                      , zorder=1.3
         )
+
+        #ax.quiver(X[::15], Z[::15], V[::15,::15], U[::15,::15], units='width',color="white")
         
     #

examples/go-cc-examples.sh

@@ -10,13 +10,13 @@ PROGRAM='scattnlay-example.bin'
 # echo Compilation done. Running...
 # time ./$PROGRAM
 
-file=test-surf-integral.cc
-echo Compile $file with gcc
-rm -f $PROGRAM
-g++ -Ofast -std=c++11 $file  ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-echo Compilation done. Running...
-./$PROGRAM
+# file=test-surf-integral.cc
+# echo Compile $file with gcc
+# rm -f $PROGRAM
+# g++ -Ofast -std=c++11 $file  ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# # g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ../src/shell-generator.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# echo Compilation done. Running...
+# ./$PROGRAM
 
 
 # file=example-minimal.cc
@@ -27,16 +27,16 @@ echo Compilation done. Running...
 # ./$PROGRAM
 
 
-# file=example-get-Mie.cc
-# echo Compile $file with gcc
-# rm -f $PROGRAM
-# # Production for multiprecision
-# #g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc read-spectra.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
+file=example-get-Mie.cc
+echo Compile $file with gcc
+rm -f $PROGRAM
+# Production for multiprecision
+#g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc read-spectra.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
 
-# # Simplified for multiprecision
-# #g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
+# Simplified for multiprecision
+#g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
 
-# # Simplified for double precision
-# g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -lrt -o $PROGRAM -march=native -mtune=native -msse4.2
-# echo Compilation done. Running...
-# ./$PROGRAM
+# Simplified for double precision
+clang++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -o $PROGRAM -march=native -mtune=native -msse4.2
+echo Compilation done. Running...
+./$PROGRAM

scattnlay.pyx

@@ -43,6 +43,17 @@ cdef inline double *npy2c(np.ndarray a):
 
 cdef extern from "py_nmie.h":
     cdef int ScattCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double anr[], double ani[], double bnr[], double bni[])
+    cdef int ExpansionCoeffs( int L,  int pl, vector[double] x,  vector[complex] m,
+                     int nmax,
+                    vector[vector[double] ]  alnr,
+                    vector[vector[double] ]  alni,
+                    vector[vector[double] ]  blnr,
+                    vector[vector[double] ]  blni,
+                    vector[vector[double] ]  clnr,
+                    vector[vector[double] ]  clni,
+                    vector[vector[double] ]  dlnr,
+                    vector[vector[double] ]  dlni)
+    # cdef int ExpansionCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double alnr[], double alni[], double blnr[], double blni[], double clnr[], double clni[], double dlnr[], double dlni[])
     cdef int nMie(int L, int pl, vector[double] x, vector[complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[])
     cdef int nField(int L, int pl, vector[double] x, vector[complex] m, int nmax, int mode_n, int mode_type, int nCoords, vector[double] Xp, vector[double] Yp, vector[double] Zp, double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[], double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[])
 
@@ -72,6 +83,45 @@ def scattcoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t
 
     return terms, an, bn
 
+def expansioncoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.int_t nmax, np.int_t pl = -1):
+    cdef Py_ssize_t i
+    cdef Py_ssize_t l
+
+    cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
+    
+    cdef np.ndarray[np.complex128_t, ndim = 3] aln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 3] bln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 3] cln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 3] dln = np.zeros((x.shape[0], x.shape[1]+1, nmax), dtype = np.complex128)
+
+    cdef np.ndarray[np.float64_t, ndim = 2] alnr
+    cdef np.ndarray[np.float64_t, ndim = 2] alni
+    cdef np.ndarray[np.float64_t, ndim = 2] blnr
+    cdef np.ndarray[np.float64_t, ndim = 2] blni
+    cdef np.ndarray[np.float64_t, ndim = 2] clnr
+    cdef np.ndarray[np.float64_t, ndim = 2] clni
+    cdef np.ndarray[np.float64_t, ndim = 2] dlnr
+    cdef np.ndarray[np.float64_t, ndim = 2] dlni
+
+    for i in range(x.shape[0]):
+        alnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        alni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        blnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        blni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        clnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        clni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        dlnr = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+        dlni = np.zeros((x.shape[1]+1,nmax), dtype = np.float64)
+
+        terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, alnr, alni, blnr, blni, clnr, clni, dlnr, dlni)
+        # terms[i] = ExpansionCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, npy2c(alnr), npy2c(alni), npy2c(blnr), npy2c(blni), npy2c(clnr), npy2c(clni), npy2c(dlnr), npy2c(dlni))
+
+        for l in range(x.shape[1]):
+            aln[l][i] = alnr[l].copy('C') + 1.0j*alni[l].copy('C')
+            bln[l][i] = blnr[l].copy('C') + 1.0j*blni[l].copy('C')
+
+    return terms, aln, bln, cln, dln
+
 def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t nmax = -1, np.int_t pl = -1):
     cdef Py_ssize_t i
 

setup.py

@@ -1,8 +1,8 @@
 #!/usr/bin/env python
 # -*- coding: UTF-8 -*-
 #
-#    Copyright (C) 2009-2017 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
-#    Copyright (C) 2013-2017 Konstantin Ladutenko <kostyfisik@gmail.com>
+#    Copyright (C) 2009-2018 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2018 Konstantin Ladutenko <kostyfisik@gmail.com>
 #
 #    This file is part of python-scattnlay
 #
@@ -26,12 +26,12 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-__version__ = '2.1'
+__version__ = '2.2' # compare with Makefile and nmie.hpp
 __title__ = 'Calculation of the scattering of EM radiation by a multilayered sphere'
 __mod__ = 'python-scattnlay'
 __author__ = 'Ovidio Peña Rodríguez'
 __email__ = 'ovidio@bytesfall.com'
-__url__ = 'http://scattering.sourceforge.net/'
+__url__ = 'https://github.com/ovidiopr/scattnlay'
 
 from distutils.core import setup
 from distutils.extension import Extension

src/nmie.cc

@@ -97,6 +97,54 @@ namespace nmie {
   }
 
   //**********************************************************************************//
+  // This function emulates a C call to calculate the scattering coefficients         //
+  // required to calculate both the near- and far-field parameters.                   //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   pl: Index of PEC layer. If there is none just send -1                          //
+  //   x: Array containing the size parameters of the layers [0..L-1]                 //
+  //   m: Array containing the relative refractive indexes of the layers [0..L-1]     //
+  //   nmax: Maximum number of multipolar expansion terms to be used for the          //
+  //         calculations. Only use it if you know what you are doing, otherwise      //
+  //         set this parameter to -1 and the function will calculate it.             //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   aln, bln ,cln, dln : Complex scattering amplitudes                             //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations in each layer   //
+  //**********************************************************************************//
+  int ExpansionCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::vector<std::complex<double> > >& aln, std::vector<std::vector<std::complex<double> > >& bln, std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
+
+    if (x.size() != L || m.size() != L)
+        throw std::invalid_argument("Declared number of layers do not fit x and m!");
+    try {
+      MultiLayerMie<FloatType> ml_mie;
+      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
+      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
+      ml_mie.SetPECLayer(pl);
+      ml_mie.SetMaxTerms(nmax);
+
+      ml_mie.calcScattCoeffs();
+      ml_mie.calcExpanCoeffs();
+
+      aln = ConvertComplexVectorVector<double>(ml_mie.GetAln());
+      bln = ConvertComplexVectorVector<double>(ml_mie.GetBln());
+      cln = ConvertComplexVectorVector<double>(ml_mie.GetCln());
+      dln = ConvertComplexVectorVector<double>(ml_mie.GetDln());
+
+      return ml_mie.GetMaxTerms();
+    } catch(const std::invalid_argument& ia) {
+      // Will catch if  ml_mie fails or other errors.
+      std::cerr << "Invalid argument: " << ia.what() << std::endl;
+      throw std::invalid_argument(ia);
+      return -1;
+    }
+    return 0;
+  }
+
+  //**********************************************************************************//
   // This function emulates a C call to calculate the actual scattering parameters    //
   // and amplitudes.                                                                  //
   //                                                                                  //

src/nmie.hpp

@@ -27,7 +27,7 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#define VERSION "2.2"  //Compare with Makefile
+#define VERSION "2.2"  //Compare with Makefile and setup.py
 #include <array>
 #include <complex>
 #include <cstdlib>
@@ -36,6 +36,7 @@
 #include <boost/math/constants/constants.hpp>
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn);
+  int ExpansionCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn, std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn);
   int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
@@ -58,6 +59,7 @@ namespace nmie {
     void RunMieCalculation();
     void RunFieldCalculation(const int mode_n=Modes::kAll, const int mode_type=Modes::kAll);
     void calcScattCoeffs();
+    void calcExpanCoeffs();
 
     // Return calculation results
     FloatType GetQext();
@@ -73,6 +75,11 @@ namespace nmie {
     std::vector<std::complex<FloatType> > GetAn(){return an_;};
     std::vector<std::complex<FloatType> > GetBn(){return bn_;};
 
+    std::vector<std::vector<std::complex<FloatType> > > GetAln(){return aln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetBln(){return bln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetCln(){return cln_;};
+    std::vector<std::vector<std::complex<FloatType> > > GetDln(){return dln_;};
+
     // Problem definition
     // Modify size of all layers
     void SetLayersSize(const std::vector<FloatType>& layer_size);
@@ -123,7 +130,6 @@ namespace nmie {
     // Scattering coefficients
     std::vector<std::complex<FloatType> > an_, bn_;
     std::vector< std::vector<std::complex<FloatType> > > aln_, bln_, cln_, dln_;
-    void calcExpanCoeffs();
     // Points for field evaluation
     std::vector< std::vector<FloatType> > coords_;
 

src/py_nmie.cc

@@ -34,17 +34,18 @@
 // Same as ScattCoeffs in 'nmie.h' but uses double arrays to return the results (useful for python).
 // This is a workaround because I have not been able to return the results using 
 // std::vector<std::complex<double> >
-int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
+int ScattCoeffs(const unsigned int L, const int pl,
+                std::vector<double>& x, std::vector<std::complex<double> >& m,
                 const int nmax, double anr[], double ani[], double bnr[], double bni[]) {
 
-  int i, result;
+  int result;
   std::vector<std::complex<double> > an, bn;
   an.resize(nmax);
   bn.resize(nmax);
 
   result = nmie::ScattCoeffs(L, pl, x, m, nmax, an, bn);
 
-  for (i = 0; i < result; i++) {
+  for (int i = 0; i < result; i++) {
     anr[i] = an[i].real();
     ani[i] = an[i].imag();
     bnr[i] = bn[i].real();
@@ -54,6 +55,52 @@ int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std:
   return result;
 }
 
+// Same as ExpansionCoeffs in 'nmie.h' but uses double arrays to
+// return the results (useful for python).  This is a workaround
+// because I have not been able to return the results using
+// std::vector<std::vector<std::complex<double> > >
+int ExpansionCoeffs(const unsigned int L, const int pl,
+                    std::vector<double>& x, std::vector<std::complex<double> >& m,
+                    const int nmax,
+                    std::vector<std::vector<double> >&  alnr,
+                    std::vector<std::vector<double> >&  alni,
+                    std::vector<std::vector<double> >&  blnr,
+                    std::vector<std::vector<double> >&  blni,
+                    std::vector<std::vector<double> >&  clnr,
+                    std::vector<std::vector<double> >&  clni,
+                    std::vector<std::vector<double> >&  dlnr,
+                    std::vector<std::vector<double> >&  dlni) {
+
+  int result;
+  std::vector< std::vector<std::complex<double> > > aln, bln, cln, dln;
+  aln.resize(L + 1);
+  bln.resize(L + 1);
+  cln.resize(L + 1);
+  dln.resize(L + 1);
+  for (unsigned int l = 0; l <= L; l++) {
+    aln[l].resize(nmax);
+    bln[l].resize(nmax);
+    cln[l].resize(nmax);
+    dln[l].resize(nmax);
+  }
+
+  result = nmie::ExpansionCoeffs(L, pl, x, m, nmax, aln, bln, cln, dln);
+
+  for (unsigned int l = 0; l <= L; l++) {
+    for (int i = 0; i < result; i++) {
+      alnr[l][i] = aln[l][i].real();
+      alni[l][i] = aln[l][i].imag();
+      blnr[l][i] = bln[l][i].real();
+      blni[l][i] = bln[l][i].imag();
+      clnr[l][i] = cln[l][i].real();
+      clni[l][i] = cln[l][i].imag();
+      dlnr[l][i] = dln[l][i].real();
+      dlni[l][i] = dln[l][i].imag();
+    }
+  }
+  return result;
+}
+
 // Same as nMie in 'nmie.h' but uses double arrays to return the results (useful for python).
 // This is a workaround because I have not been able to return the results using 
 // std::vector<std::complex<double> >

src/py_nmie.h

@@ -28,8 +28,27 @@
 #include <complex>
 #include <vector>
 
-int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
-                const int nmax, double anr[], double ani[], double bnr[], double bni[]);
+int ScattCoeffs(const unsigned int L, const int pl,
+                std::vector<double>& x, std::vector<std::complex<double> >& m,
+                const int nmax,
+                double anr[], double ani[], double bnr[], double bni[]);
+
+int ExpansionCoeffs(const unsigned int L, const int pl,
+                    std::vector<double>& x, std::vector<std::complex<double> >& m,
+                    const int nmax,
+                    std::vector<std::vector<double> >&  alnr,
+                    std::vector<std::vector<double> >&  alni,
+                    std::vector<std::vector<double> >&  blnr,
+                    std::vector<std::vector<double> >&  blni,
+                    std::vector<std::vector<double> >&  clnr,
+                    std::vector<std::vector<double> >&  clni,
+                    std::vector<std::vector<double> >&  dlnr,
+                    std::vector<std::vector<double> >&  dlni);
+/* int ExpansionCoeffs(const unsigned int L, const int pl, */
+/*                     std::vector<double>& x, std::vector<std::complex<double> >& m, */
+/*                     const int nmax, */
+/*                     double alnr[], double alni[], double blnr[], double blni[], */
+/*                     double clnr[], double clni[], double dlnr[], double dlni[]); */
 
 int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
          const int nTheta, std::vector<double>& Theta, const int nmax,