Several changes to complete the features

farfield.cc

@@ -69,29 +69,23 @@ int main(int argc, char *argv[]) {
     // for (auto arg : args) std::cout<< arg <<std::endl;
     std::string comment;
     int has_comment = 0;
-    int i;
     unsigned int L = 0;
     std::vector<double> x, Theta;
     std::vector<std::complex<double> > m, S1, S2;
     double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo;
 
-    std::vector<std::complex<double> > mw, S1w, S2w;
-    double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow;
-
     double ti = 0.0, tf = 90.0;
-    int nt = 0;    
+    int nt = 0;
     if (argc < 5) throw std::invalid_argument(error_msg);
-    
-    //strcpy(comment, "");
-    // for (i = 1; i < argc; i++) {
-    int mode = -1; 
+
+    int mode = -1;
     double tmp_mr;
     for (auto arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
 
-      // Detecting new read mode (if it is a valid -key) 
+      // Detecting new read mode (if it is a valid -key)
       if (arg == "-l") {
         mode = read_L;
         continue;
@@ -154,8 +148,6 @@ int main(int argc, char *argv[]) {
         Theta.resize(nt);
         S1.resize(nt);
         S2.resize(nt);
-        S1w.resize(nt);
-        S2w.resize(nt);
         continue;
       }
 
@@ -166,35 +158,34 @@ int main(int argc, char *argv[]) {
       }
     }
 
-    if ( (x.size() != m.size()) || (L != x.size()) ) 
+    if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n") + error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) ) 
+    if ( (0 == m.size()) || ( 0 == x.size()) )
       throw std::invalid_argument(std::string("Empty structure!\n") + error_msg);
-    
+
     if (nt < 0) {
       printf("Error reading Theta.\n");
       return -1;
     } else if (nt == 1) {
       Theta[0] = ti*PI/180.0;
     } else {
-      for (i = 0; i < nt; i++) {
-      Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
+      for (int i = 0; i < nt; i++) {
+        Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
       }
     }
 
-    nmie::nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-   
+    nmie::nMie(L, -1, x, m, nt, Theta, -1, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
 
     if (has_comment) {
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", comment.c_str(), Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
     } else {
       printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
     }
-    
+
     if (nt > 0) {
       printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
-      
-      for (i = 0; i < nt; i++) {
+
+      for (int i = 0; i < nt; i++) {
         printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e\n", Theta[i]*180.0/PI, S1[i].real(), S1[i].imag(), S2[i].real(), S2[i].imag());
       }
     }
@@ -204,7 +195,7 @@ int main(int argc, char *argv[]) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
-  }  
+  }
     return 0;
 }
 

nearfield.cc

@@ -0,0 +1,265 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program 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.                                           //
+//                                                                                  //
+//    This program 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.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//                                                                                  //
+//    You should have received a copy of the GNU General Public License             //
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+
+#include <algorithm>
+#include <complex>
+#include <functional>
+#include <iostream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+#include <string.h>
+#include "nmie.h"
+
+const double PI=3.14159265358979323846;
+
+//***********************************************************************************//
+// This is the main function of 'scattnlay', here we read the parameters as          //
+// arguments passed to the program which should be executed with the following       //
+// syntaxis:                                                                         //
+// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]  //
+//                                                                                   //
+// When all the parameters were correctly passed we setup the integer L (the         //
+// number of layers) and the arrays x and m, containing the size parameters and      //
+// refractive indexes of the layers, respectively and call the function nMie.        //
+// If the calculation is successful the results are printed with the following       //
+// format:                                                                           //
+//                                                                                   //
+//    * If no comment was passed:                                                    //
+//        'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                                    //
+//                                                                                   //
+//    * If a comment was passed:                                                     //
+//        'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                           //
+//***********************************************************************************//
+int main(int argc, char *argv[]) {
+  try {
+    std::vector<std::string> args;
+    args.assign(argv, argv + argc);
+    std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + args[0]
+                          + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] "
+                          + "[-p xi xf nx yi yf ny zi zf nz] [-c comment]\n");
+    enum mode_states {read_L, read_x, read_mr, read_mi, read_xi, read_xf, read_nx, read_yi, read_yf, read_ny, read_zi, read_zf, read_nz, read_comment};
+    // for (auto arg : args) std::cout<< arg <<std::endl;
+    std::string comment;
+    int has_comment = 0;
+    unsigned int L = 0;
+    std::vector<double> x, Xp, Yp, Zp;
+    std::vector<std::complex<double> > m;
+    std::vector<std::vector<std::complex<double> > > E, H;
+
+    double xi = 0.0, xf = 0.0, yi = 0.0, yf = 0.0, zi = 0.0, zf = 0.0;
+    double dx = 0.0, dy = 0.0, dz = 0.0;
+    int nx = 0, ny = 0, nz = 0;
+    long total_points = 0;
+    if (argc < 5) throw std::invalid_argument(error_msg);
+
+    int mode = -1;
+    double tmp_mr;
+    for (auto arg : args) {
+      // For each arg in args list we detect the change of the current
+      // read mode or read the arg. The reading args algorithm works
+      // as a finite-state machine.
+
+      // Detecting new read mode (if it is a valid -key)
+      if (arg == "-l") {
+        mode = read_L;
+        continue;
+      }
+
+      if (arg == "-p") {
+        if ((mode != read_x) && (mode != read_comment))
+          throw std::invalid_argument(std::string("Unfinished layer!\n") + error_msg);
+        mode = read_xi;
+        continue;
+      }
+
+      if (arg == "-c") {
+        if ((mode != read_x) && (mode != read_nz))
+          throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg);
+        mode = read_comment;
+        continue;
+      }
+
+      // Reading data. For invalid date the exception will be thrown
+      // with the std:: and catched in the end.
+      if (mode == read_L) {
+        L = std::stoi(arg);
+        mode = read_x;
+        continue;
+      }
+
+      if (mode == read_x) {
+        x.push_back(std::stod(arg));
+        mode = read_mr;
+        continue;
+      }
+
+      if (mode == read_mr) {
+        tmp_mr = std::stod(arg);
+        mode = read_mi;
+        continue;
+      }
+
+      if (mode == read_mi) {
+        m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
+        mode = read_x;
+        continue;
+      }
+
+      if (mode == read_xi) {
+        xi = std::stod(arg);
+        mode = read_xf;
+        continue;
+      }
+
+      if (mode == read_xf) {
+        xf = std::stod(arg);
+        mode = read_nx;
+        continue;
+      }
+
+      if (mode == read_nx) {
+        nx = std::stoi(arg);
+        mode = read_yi;
+        continue;
+      }
+
+      if (mode == read_yi) {
+        yi = std::stod(arg);
+        mode = read_yf;
+        continue;
+      }
+
+      if (mode == read_yf) {
+        yf = std::stod(arg);
+        mode = read_ny;
+        continue;
+      }
+
+      if (mode == read_ny) {
+        ny = std::stoi(arg);
+        mode = read_zi;
+        continue;
+      }
+
+      if (mode == read_zi) {
+        zi = std::stod(arg);
+        mode = read_zf;
+        continue;
+      }
+
+      if (mode == read_zf) {
+        zf = std::stod(arg);
+        mode = read_nz;
+        continue;
+      }
+
+      if (mode == read_nz) {
+        nz = std::stoi(arg);
+        total_points = nx*ny*nz;
+        if (total_points <= 0)
+          throw std::invalid_argument(std::string("Nothing to do! You must define the grid to calculate the fields.\n") + error_msg);
+
+        Xp.resize(total_points);
+        Yp.resize(total_points);
+        Zp.resize(total_points);
+
+        E.resize(total_points);
+        H.resize(total_points);
+        for (long i = 0; i < total_points; i++) {
+          E[i].resize(3);
+          H[i].resize(3);
+        }
+        continue;
+      }
+
+      if (mode ==  read_comment) {
+        comment = arg;
+        has_comment = 1;
+        continue;
+      }
+    }
+
+    if ( (x.size() != m.size()) || (L != x.size()) )
+      throw std::invalid_argument(std::string("Broken structure!\n") + error_msg);
+    if ( (0 == m.size()) || ( 0 == x.size()) )
+      throw std::invalid_argument(std::string("Empty structure!\n") + error_msg);
+
+    if (nx == 1)
+      dx = 0.0;
+    else
+      dx = (xf - xi)/(nx - 1);
+
+    if (ny == 1)
+      dy = 0.0;
+    else
+      dy = (yf - yi)/(ny - 1);
+
+    if (nz == 1)
+      dz = 0.0;
+    else
+      dz = (zf - zi)/(nz - 1);
+
+    for (int i = 0; i < nx; i++) {
+      for (int j = 0; j < ny; j++) {
+        for (int k = 0; k < nz; k++) {
+          Xp[i*ny + j*nz + k] = xi + (double)i*dx;
+          Yp[i*ny + j*nz + k] = yi + (double)j*dy;
+          Zp[i*ny + j*nz + k] = zi + (double)k*dz;
+        }
+      }
+    }
+
+    nmie::nField(L, -1, x, m, -1, total_points, Xp, Yp, Zp, E, H);
+
+    if (has_comment)
+      printf("%6s\n", comment.c_str());
+
+    if (total_points > 0) {
+      printf("         X,          Y,          Z,         Ex.r,         Ex.i,         Ey.r,         Ey.i,         Ez.r,         Ez.i,         Hx.r,         Hx.i,         Hy.r,         Hy.i,         Hz.r,         Hz.i\n");
+
+      for (long i = 0; i < total_points; i++) {
+        printf("%10.7f, %10.7f, %10.7f, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n",
+               Xp[i], Yp[i], Zp[i],
+               E[i][0].real(), E[i][0].imag(), E[i][1].real(), E[i][1].imag(), E[i][2].real(), E[i][2].imag(),
+               H[i][0].real(), H[i][0].imag(), H[i][1].real(), H[i][1].imag(), H[i][2].real(), H[i][2].imag());
+      }
+    }
+
+
+  } catch( const std::invalid_argument& ia ) {
+    // Will catch if  multi_layer_mie fails or other errors.
+    std::cerr << "Invalid argument: " << ia.what() << std::endl;
+    return -1;
+  }
+    return 0;
+}
+
+

nmie.cc

@@ -54,6 +54,51 @@ 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:                                                               //
+  //   an, bn: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  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) {
+
+    if (x.size() != L || m.size() != L)
+        throw std::invalid_argument("Declared number of layers do not fit x and m!");
+    try {
+      MultiLayerMie ml_mie;
+      ml_mie.SetLayersSize(x);
+      ml_mie.SetLayersIndex(m);
+      ml_mie.SetPECLayer(pl);
+      ml_mie.SetMaxTerms(nmax);
+
+      ml_mie.calcScattCoeffs();
+
+      an = ml_mie.GetAn();
+      bn = ml_mie.GetBn();
+
+      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.                                                                  //
   //                                                                                  //
@@ -89,26 +134,28 @@ namespace nmie {
     if (Theta.size() != nTheta)
         throw std::invalid_argument("Declared number of sample for Theta is not correct!");
     try {
-      MultiLayerMie multi_layer_mie;
-      multi_layer_mie.SetLayersSize(x);
-      multi_layer_mie.SetLayersIndex(m);
-      multi_layer_mie.SetAngles(Theta);
-      multi_layer_mie.SetPECLayer(pl);
-      multi_layer_mie.SetMaxTerms(nmax);
-
-      multi_layer_mie.RunMieCalculation();
-
-      *Qext = multi_layer_mie.GetQext();
-      *Qsca = multi_layer_mie.GetQsca();
-      *Qabs = multi_layer_mie.GetQabs();
-      *Qbk = multi_layer_mie.GetQbk();
-      *Qpr = multi_layer_mie.GetQpr();
-      *g = multi_layer_mie.GetAsymmetryFactor();
-      *Albedo = multi_layer_mie.GetAlbedo();
-      S1 = multi_layer_mie.GetS1();
-      S2 = multi_layer_mie.GetS2();
+      MultiLayerMie ml_mie;
+      ml_mie.SetLayersSize(x);
+      ml_mie.SetLayersIndex(m);
+      ml_mie.SetAngles(Theta);
+      ml_mie.SetPECLayer(pl);
+      ml_mie.SetMaxTerms(nmax);
+
+      ml_mie.RunMieCalculation();
+
+      *Qext = ml_mie.GetQext();
+      *Qsca = ml_mie.GetQsca();
+      *Qabs = ml_mie.GetQabs();
+      *Qbk = ml_mie.GetQbk();
+      *Qpr = ml_mie.GetQpr();
+      *g = ml_mie.GetAsymmetryFactor();
+      *Albedo = ml_mie.GetAlbedo();
+      S1 = ml_mie.GetS1();
+      S2 = ml_mie.GetS2();
+
+      return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {
-      // Will catch if  multi_layer_mie fails or other errors.
+      // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
       return -1;
@@ -250,16 +297,18 @@ namespace nmie {
       if (f.size() != 3)
         throw std::invalid_argument("Field H is not 3D!");
     try {
-      MultiLayerMie multi_layer_mie;
-      //multi_layer_mie.SetPECLayer(pl); // TODO add PEC layer to field plotting
-      multi_layer_mie.SetLayersSize(x);
-      multi_layer_mie.SetLayersIndex(m);
-      multi_layer_mie.SetFieldCoords({Xp_vec, Yp_vec, Zp_vec});
-      multi_layer_mie.RunFieldCalculation();
-      E = multi_layer_mie.GetFieldE();
-      H = multi_layer_mie.GetFieldH();
+      MultiLayerMie ml_mie;
+      //ml_mie.SetPECLayer(pl); // TODO add PEC layer to field plotting
+      ml_mie.SetLayersSize(x);
+      ml_mie.SetLayersIndex(m);
+      ml_mie.SetFieldCoords({Xp_vec, Yp_vec, Zp_vec});
+      ml_mie.RunFieldCalculation();
+      E = ml_mie.GetFieldE();
+      H = ml_mie.GetFieldH();
+
+      return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {
-      // Will catch if  multi_layer_mie fails or other errors.
+      // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
       return - 1;
@@ -723,7 +772,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  void MultiLayerMie::ScattCoeffs() {
+  void MultiLayerMie::calcScattCoeffs() {
 
     isScaCoeffsCalc_ = false;
 
@@ -873,7 +922,7 @@ namespace nmie {
       }
     }  // end of for an and bn terms
     isScaCoeffsCalc_ = true;
-  }  // end of MultiLayerMie::ScattCoeffs(...)
+  }  // end of MultiLayerMie::calcScattCoeffs()
 
 
   //**********************************************************************************//
@@ -917,7 +966,7 @@ namespace nmie {
     isMieCalculated_ = false;
 
     // Calculate scattering coefficients
-    ScattCoeffs();
+    calcScattCoeffs();
 
     if (!isScaCoeffsCalc_) // TODO seems to be unreachable
       throw std::invalid_argument("Calculation of scattering coefficients failed!");
@@ -996,7 +1045,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  void MultiLayerMie::ExpanCoeffs() {
+  void MultiLayerMie::calcExpanCoeffs() {
     if (!isScaCoeffsCalc_)
       throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
 
@@ -1086,7 +1135,7 @@ namespace nmie {
     }
 
     isExpCoeffsCalc_ = true;
-  }  // end of   void MultiLayerMie::ExpanCoeffs()
+  }  // end of   void MultiLayerMie::calcExpanCoeffs()
 
 
   //**********************************************************************************//
@@ -1194,10 +1243,10 @@ namespace nmie {
     double Rho, Theta, Phi;
 
     // Calculate scattering coefficients an_ and bn_
-    ScattCoeffs();
+    calcScattCoeffs();
 
     // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
-    ExpanCoeffs();
+    calcExpanCoeffs();
 
     long total_points = coords_[0].size();
     E_.resize(total_points);

nmie.h

@@ -33,7 +33,7 @@
 #include <vector>
 
 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 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 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);
@@ -45,6 +45,7 @@ namespace nmie {
     // Run calculation
     void RunMieCalculation();
     void RunFieldCalculation();
+    void calcScattCoeffs();
 
     // Return calculation results
     double GetQext();
@@ -122,8 +123,7 @@ namespace nmie {
                        const double& Pi, const double& Tau, const double& n,
                        std::vector<std::complex<double> >& Mo1n, std::vector<std::complex<double> >& Me1n, 
                        std::vector<std::complex<double> >& No1n, std::vector<std::complex<double> >& Ne1n);
-    void ScattCoeffs();
-    void ExpanCoeffs();
+    void calcExpanCoeffs();
 
     void calcField(const double Rho, const double Theta, const double Phi,
                    std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);

py_nmie.cc

@@ -31,6 +31,29 @@
 #include "nmie.h"
 #include "py_nmie.h"
 
+// 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,
+                const int nmax, double anr[], double ani[], double bnr[], double bni[]) {
+
+  int i, 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++) {
+    anr[i] = an[i].real();
+    ani[i] = an[i].imag();
+    bnr[i] = bn[i].real();
+    bni[i] = bn[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> >

py_nmie.h

@@ -28,6 +28,9 @@
 #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 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,
          double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,

scattnlay.pyx

@@ -42,10 +42,37 @@ cdef inline double *npy2c(np.ndarray a):
     return <double *>(a.data)
 
 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 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 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[])
 
-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 pl = -1, np.int_t nmax = -1):
+def scattcoeffs(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 np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
+
+    cdef np.ndarray[np.complex128_t, ndim = 2] an = np.zeros((x.shape[0], nmax), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 2] bn = np.zeros((x.shape[0], nmax), dtype = np.complex128)
+
+    cdef np.ndarray[np.float64_t, ndim = 1] anr
+    cdef np.ndarray[np.float64_t, ndim = 1] ani
+    cdef np.ndarray[np.float64_t, ndim = 1] bnr
+    cdef np.ndarray[np.float64_t, ndim = 1] bni
+
+    for i in range(x.shape[0]):
+        anr = np.zeros(nmax, dtype = np.float64)
+        ani = np.zeros(nmax, dtype = np.float64)
+        bnr = np.zeros(nmax, dtype = np.float64)
+        bni = np.zeros(nmax, dtype = np.float64)
+
+        terms[i] = ScattCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, npy2c(anr), npy2c(ani), npy2c(bnr), npy2c(bni))
+
+        an[i] = anr.copy('C') + 1.0j*ani.copy('C')
+        bn[i] = bnr.copy('C') + 1.0j*bni.copy('C')
+
+    return terms, an, bn
+
+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
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
@@ -79,8 +106,7 @@ def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t,
 
     return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
 
-#def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords = np.zeros((0, 3), dtype = np.float64), np.int_t pl = 0, np.int_t nmax = 0):
-def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords, np.int_t pl = 0, np.int_t nmax = 0):
+def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords, np.int_t nmax = -1, np.int_t pl = -1):
     cdef Py_ssize_t i
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)

setup_cython.py

@@ -4,7 +4,7 @@
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
 #    Copyright (C) 2013-2015 Konstantin Ladutenko <kostyfisik@gmail.com>
 #
-#    This file is part of python-scattnlay
+#    This file is part of scattnlay
 #
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by

tests/python/field.py → tests/python/field-nanoshell.py

@@ -2,8 +2,9 @@
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015 Konstantin Ladutenko <kostyfisik@gmail.com>
 #
-#    This file is part of python-scattnlay
+#    This file is part of scattnlay
 #
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -26,18 +27,11 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # This test case calculates the electric field in the 
-# XY plane, for an spherical silver nanoparticle
-# embedded in glass.
+# XY plane, for a silver nanoshell embedded in water.
 
-# Refractive index values correspond to a wavelength of
-# 400 nm. Maximum of the surface plasmon resonance (and,
-# hence, of electric field) is expected under those
-# conditions.
-import scattnlay
-
-#import os
-#path = os.path.dirname(scattnlay.__file__)
-#print(scattnlay.__file__)
+# Refractive index values correspond to the wavelength
+# where maximum of the surface plasmon resonance (and,
+# hence, of electric field) is expected.
 
 from scattnlay import fieldnlay
 import numpy as np
@@ -132,7 +126,7 @@ try:
     plt.clf()
     plt.close()
 finally:
-    np.savetxt("field.txt", result, fmt = "%.5f")
+    np.savetxt("field-nanoshell.txt", result, fmt = "%.5f")
     print result
 
 

tests/python/test01.py

@@ -50,23 +50,20 @@
 import scattnlay
 
 import os
-path = os.path.dirname(scattnlay.__file__)
-print(scattnlay.__file__)
 
 from scattnlay import scattnlay
 import numpy as np
-# import os
-# import inspect
-# inspect.getfile(scattnlay)
-
-x = np.ones((400, 5), dtype = np.float64)
-x[:, 4] = np.arange(0.25, 100.25, 0.25)
-x[:, 0] = 0.1**(1.0/3.0)*x[:, 4]
-x[:, 1] = 0.36**(1.0/3.0)*x[:, 4]
-x[:, 2] = 0.404**(1.0/3.0)*x[:, 4]
-x[:, 3] = 0.7706**(1.0/3.0)*x[:, 4]
-
-m = np.ones((400, 5), dtype = np.complex128)
+
+size = np.arange(0.25, 100.25, 0.25)
+
+x = np.ones((len(size), 5), dtype = np.float64)
+x[:, 0] = 0.1**(1.0/3.0)*size
+x[:, 1] = 0.36**(1.0/3.0)*size
+x[:, 2] = 0.404**(1.0/3.0)*size
+x[:, 3] = 0.7706**(1.0/3.0)*size
+x[:, 4] = size
+
+m = np.ones((len(size), 5), dtype = np.complex128)
 m[:, 0] *= 1.8 + 1.7j
 m[:, 1] *= 0.8 + 0.7j
 m[:, 2] *= 1.2 + 0.09j
@@ -95,7 +92,7 @@ try:
 
     plt.xlabel('X')
     
-    #plt.show()
+    plt.show()
 finally:
     np.savetxt("test01.txt", result, fmt = "%.5f")
     print result

tests/shell/field-nanoshell.sh

@@ -0,0 +1,37 @@
+#! /bin/sh
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015 Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of scattnlay
+#
+#    This program 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.
+#
+#    This program 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.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite the following reference:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# This test case calculates the electric field in the 
+# XY plane, for a silver nanoshell embedded in water.
+
+# Refractive index values correspond to the wavelength
+# where maximum of the surface plasmon resonance (and,
+# hence, of electric field) is expected.
+
+PROGRAM='../../../fieldnlay'
+
+$PROGRAM -l 2 0.38989409 1.16177963 0.00000000 0.46787291 0.42850284 5.47718289 -p -1.55957635 1.55957635 501 -1.55957635 1.55957635 501 0.00000000 0.00000000 1