Merge branch 'master' of github.com:ovidiopr/scattnlay

Makefile

@@ -7,19 +7,26 @@ VERSION=0.3.1
 all:
 	@echo "make source - Create source package"
 	@echo "make cython - Convert Cython code to c++"
+	@echo "make python_ext - Create Python extension using C++ code"
+	@echo "make cython_ext - Create Python extension using Cython code"
 	@echo "make install - Install on local system"
 	@echo "make buildrpm - Generate a rpm package"
 	@echo "make builddeb - Generate a deb package"
 	@echo "make standalone - Create a standalone program"
 	@echo "make clean - Delete temporal files"
-	make standalone
+#	make standalone
 
 source:
 	$(PYTHON) setup.py sdist $(COMPILE) --dist-dir=../
 
 cython: scattnlay.pyx
 	cython --cplus scattnlay.pyx
-	mv scattnlay.cpp scattnlay.cc
+
+python_ext: nmie.cc py_nmie.cc scattnlay.cpp
+	export CFLAGS='-std=c++11' && python setup.py build_ext --inplace
+
+cython_ext: nmie.cc py_nmie.cc scattnlay.pyx
+	export CFLAGS='-std=c++11' && python setup_cython.py build_ext --inplace
 
 install:
 	$(PYTHON) setup.py install --root $(DESTDIR) $(COMPILE)
@@ -37,7 +44,7 @@ builddeb:
 	dpkg-buildpackage -i -I -rfakeroot
 
 standalone: standalone.cc nmie.cc
-	c++ -DNDEBUG -O2 -std=c++11 standalone.cc nmie.cc nmie-wrapper.cc -lm -o scattnlay
+	c++ -DNDEBUG -O2 -std=c++11 standalone.cc nmie.cc -lm -o scattnlay
 	mv scattnlay ../
 
 clean:

go.sh

@@ -10,7 +10,7 @@ rm -f ../scattnlay
 
 #google profiler  ######## FAST!!!
 echo Uncomment next line to compile compare.cc
-#g++ -Ofast -std=c++11 $file nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.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 nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.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 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay-g.bin -ltcmalloc -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -g
 

nmie-old.cc

@@ -0,0 +1,1005 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.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 library implements the algorithm for a multilayered sphere described by:    //
+//    [1] W. Yang, "Improved recursive algorithm for light scattering by a          //
+//        multilayered sphere,” Applied Optics,  vol. 42, Mar. 2003, pp. 1710-1720. //
+//                                                                                  //
+// You can find the description of all the used equations in:                       //
+//    [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//                                                                                  //
+// Hereinafter all equations numbers refer to [2]                                   //
+//**********************************************************************************//
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "nmie.h"
+
+#define round(x) ((x) >= 0 ? (int)((x) + 0.5):(int)((x) - 0.5))
+
+const double PI=3.14159265358979323846;
+// light speed [m s-1]
+double const cc = 2.99792458e8;
+// assume non-magnetic (MU=MU0=const) [N A-2]
+double const mu = 4.0*PI*1.0e-7;
+
+// Calculate Nstop - equation (17)
+int Nstop(double xL) {
+  int result;
+
+  if (xL <= 8) {
+    result = round(xL + 4*pow(xL, 1.0/3.0) + 1);
+  } else if (xL <= 4200) {
+    result = round(xL + 4.05*pow(xL, 1.0/3.0) + 2);
+  } else {
+    result = round(xL + 4*pow(xL, 1.0/3.0) + 2);
+  }
+
+  return result;
+}
+
+//**********************************************************************************//
+int Nmax(int L, int fl, int pl,
+         std::vector<double> x,
+         std::vector<std::complex<double> > m) {
+  int i, result, ri, riM1;
+  result = Nstop(x[L - 1]);
+  for (i = fl; i < L; i++) {
+    if (i > pl) {
+      ri = round(std::abs(x[i]*m[i]));
+    } else {
+      ri = 0;
+    }
+    if (result < ri) {
+      result = ri;
+    }
+
+    if ((i > fl) && ((i - 1) > pl)) {
+      riM1 = round(std::abs(x[i - 1]* m[i]));
+    } else {
+      riM1 = 0;
+    }
+    if (result < riM1) {
+      result = riM1;
+    }
+  }
+  return result + 15;
+}
+
+//**********************************************************************************//
+// This function calculates the spherical Bessel (jn) and Hankel (h1n) functions    //
+// and their derivatives for a given complex value z. See pag. 87 B&H.              //
+//                                                                                  //
+// Input parameters:                                                                //
+//   z: Real argument to evaluate jn and h1n                                        //
+//   nmax: Maximum number of terms to calculate jn and h1n                          //
+//                                                                                  //
+// Output parameters:                                                               //
+//   jn, h1n: Spherical Bessel and Hankel functions                                 //
+//   jnp, h1np: Derivatives of the spherical Bessel and Hankel functions            //
+//                                                                                  //
+// The implementation follows the algorithm by I.J. Thompson and A.R. Barnett,      //
+// Comp. Phys. Comm. 47 (1987) 245-257.                                             //
+//                                                                                  //
+// Complex spherical Bessel functions from n=0..nmax-1 for z in the upper half      //
+// plane (Im(z) > -3).                                                              //
+//                                                                                  //
+//     j[n]   = j/n(z)                Regular solution: j[0]=sin(z)/z               //
+//     j'[n]  = d[j/n(z)]/dz                                                        //
+//     h1[n]  = h[0]/n(z)             Irregular Hankel function:                    //
+//     h1'[n] = d[h[0]/n(z)]/dz                h1[0] = j0(z) + i*y0(z)              //
+//                                                   = (sin(z)-i*cos(z))/z          //
+//                                                   = -i*exp(i*z)/z                //
+// Using complex CF1, and trigonometric forms for n=0 solutions.                    //
+//**********************************************************************************//
+int sbesjh(std::complex<double> z, int nmax, std::vector<std::complex<double> >& jn, std::vector<std::complex<double> >& jnp, std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np) {
+
+  const int limit = 20000;
+  double const accur = 1.0e-12;
+  double const tm30 = 1e-30;
+
+  int n;
+  double absc;
+  std::complex<double> zi, w;
+  std::complex<double> pl, f, b, d, c, del, jn0, jndb, h1nldb, h1nbdb;
+
+  absc = std::abs(std::real(z)) + std::abs(std::imag(z));
+  if ((absc < accur) || (std::imag(z) < -3.0)) {
+    return -1;
+  }
+
+  zi = 1.0/z;
+  w = zi + zi;
+
+  pl = double(nmax)*zi;
+
+  f = pl + zi;
+  b = f + f + zi;
+  d = 0.0;
+  c = f;
+  for (n = 0; n < limit; n++) {
+    d = b - d;
+    c = b - 1.0/c;
+
+    absc = std::abs(std::real(d)) + std::abs(std::imag(d));
+    if (absc < tm30) {
+      d = tm30;
+    }
+
+    absc = std::abs(std::real(c)) + std::abs(std::imag(c));
+    if (absc < tm30) {
+      c = tm30;
+    }
+
+    d = 1.0/d;
+    del = d*c;
+    f = f*del;
+    b += w;
+
+    absc = std::abs(std::real(del - 1.0)) + std::abs(std::imag(del - 1.0));
+
+    if (absc < accur) {
+      // We have obtained the desired accuracy
+      break;
+    }
+  }
+
+  if (absc > accur) {
+    // We were not able to obtain the desired accuracy
+    return -2;
+  }
+
+  jn[nmax - 1] = tm30;
+  jnp[nmax - 1] = f*jn[nmax - 1];
+
+  // Downward recursion to n=0 (N.B.  Coulomb Functions)
+  for (n = nmax - 2; n >= 0; n--) {
+    jn[n] = pl*jn[n + 1] + jnp[n + 1];
+    jnp[n] = pl*jn[n] - jn[n + 1];
+    pl = pl - zi;
+  }
+
+  // Calculate the n=0 Bessel Functions
+  jn0 = zi*std::sin(z);
+  h1n[0] = std::exp(std::complex<double>(0.0, 1.0)*z)*zi*(-std::complex<double>(0.0, 1.0));
+  h1np[0] = h1n[0]*(std::complex<double>(0.0, 1.0) - zi);
+
+  // Rescale j[n], j'[n], converting to spherical Bessel functions.
+  // Recur   h1[n], h1'[n] as spherical Bessel functions.
+  w = 1.0/jn[0];
+  pl = zi;
+  for (n = 0; n < nmax; n++) {
+    jn[n] = jn0*(w*jn[n]);
+    jnp[n] = jn0*(w*jnp[n]) - zi*jn[n];
+    if (n != 0) {
+      h1n[n] = (pl - zi)*h1n[n - 1] - h1np[n - 1];
+
+      // check if hankel is increasing (upward stable)
+      if (std::abs(h1n[n]) < std::abs(h1n[n - 1])) {
+        jndb = z;
+        h1nldb = h1n[n];
+        h1nbdb = h1n[n - 1];
+      }
+
+      pl += zi;
+
+      h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
+    }
+  }
+
+  // success
+  return 0;
+}
+
+//**********************************************************************************//
+// This function calculates the spherical Bessel functions (bj and by) and the      //
+// logarithmic derivative (bd) for a given complex value z. See pag. 87 B&H.        //
+//                                                                                  //
+// Input parameters:                                                                //
+//   z: Complex argument to evaluate bj, by and bd                                  //
+//   nmax: Maximum number of terms to calculate bj, by and bd                       //
+//                                                                                  //
+// Output parameters:                                                               //
+//   bj, by: Spherical Bessel functions                                             //
+//   bd: Logarithmic derivative                                                     //
+//**********************************************************************************//
+void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<double> >& bj, std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd) {
+
+    std::vector<std::complex<double> > jn, jnp, h1n, h1np;
+    jn.resize(nmax);
+    jnp.resize(nmax);
+    h1n.resize(nmax);
+    h1np.resize(nmax);
+
+    // TODO verify that the function succeeds
+    int ifail = sbesjh(z, nmax, jn, jnp, h1n, h1np);
+
+    for (int n = 0; n < nmax; n++) {
+      bj[n] = jn[n];
+      by[n] = (h1n[n] - jn[n])/std::complex<double>(0.0, 1.0);
+      bd[n] = jnp[n]/jn[n] + 1.0/z;
+    }
+}
+
+// external scattering field = incident + scattered
+// BH p.92 (4.37), 94 (4.45), 95 (4.50)
+// assume: medium is non-absorbing; refim = 0; Uabs = 0
+void fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double> Pi, std::vector<double> Tau,
+              std::vector<std::complex<double> > an, std::vector<std::complex<double> > bn,
+              std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
+
+  int i, n, n1;
+  double rn;
+  std::complex<double> ci, zn, xxip, encap;
+  std::vector<std::complex<double> > vm3o1n, vm3e1n, vn3o1n, vn3e1n;
+  vm3o1n.resize(3);
+  vm3e1n.resize(3);
+  vn3o1n.resize(3);
+  vn3e1n.resize(3);
+
+  std::vector<std::complex<double> > Ei, Hi, Es, Hs;
+  Ei.resize(3);
+  Hi.resize(3);
+  Es.resize(3);
+  Hs.resize(3);
+  for (i = 0; i < 3; i++) {
+    Ei[i] = std::complex<double>(0.0, 0.0);
+    Hi[i] = std::complex<double>(0.0, 0.0);
+    Es[i] = std::complex<double>(0.0, 0.0);
+    Hs[i] = std::complex<double>(0.0, 0.0);
+  }
+
+  std::vector<std::complex<double> > bj, by, bd;
+  bj.resize(nmax+1);
+  by.resize(nmax+1);
+  bd.resize(nmax+1);
+
+  // Calculate spherical Bessel and Hankel functions
+  sphericalBessel(Rho, nmax, bj, by, bd);
+
+  ci = std::complex<double>(0.0, 1.0);
+  for (n = 0; n < nmax; n++) {
+    n1 = n + 1;
+    rn = double(n + 1);
+
+    zn = bj[n1] + ci*by[n1];
+    xxip = Rho*(bj[n] + ci*by[n]) - rn*zn;
+
+    vm3o1n[0] = std::complex<double>(0.0, 0.0);
+    vm3o1n[1] = std::cos(Phi)*Pi[n]*zn;
+    vm3o1n[2] = -std::sin(Phi)*Tau[n]*zn;
+    vm3e1n[0] = std::complex<double>(0.0, 0.0);
+    vm3e1n[1] = -std::sin(Phi)*Pi[n]*zn;
+    vm3e1n[2] = -std::cos(Phi)*Tau[n]*zn;
+    vn3o1n[0] = std::sin(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
+    vn3o1n[1] = std::sin(Phi)*Tau[n]*xxip/Rho;
+    vn3o1n[2] = std::cos(Phi)*Pi[n]*xxip/Rho;
+    vn3e1n[0] = std::cos(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
+    vn3e1n[1] = std::cos(Phi)*Tau[n]*xxip/Rho;
+    vn3e1n[2] = -std::sin(Phi)*Pi[n]*xxip/Rho;
+
+    // scattered field: BH p.94 (4.45)
+    encap = std::pow(ci, rn)*(2.0*rn + 1.0)/(rn*rn + rn);
+    for (i = 0; i < 3; i++) {
+      Es[i] = Es[i] + encap*(ci*an[n]*vn3e1n[i] - bn[n]*vm3o1n[i]);
+      Hs[i] = Hs[i] + encap*(ci*bn[n]*vn3o1n[i] + an[n]*vm3e1n[i]);
+    }
+  }
+
+  // incident E field: BH p.89 (4.21); cf. p.92 (4.37), p.93 (4.38)
+  // basis unit vectors = er, etheta, ephi
+  std::complex<double> eifac = std::exp(std::complex<double>(0.0, Rho*std::cos(Theta)));
+
+  Ei[0] = eifac*std::sin(Theta)*std::cos(Phi);
+  Ei[1] = eifac*std::cos(Theta)*std::cos(Phi);
+  Ei[2] = -eifac*std::sin(Phi);
+
+  // magnetic field
+  double hffact = 1.0/(cc*mu);
+  for (i = 0; i < 3; i++) {
+    Hs[i] = hffact*Hs[i];
+  }
+
+  // incident H field: BH p.26 (2.43), p.89 (4.21)
+  std::complex<double> hffacta = hffact;
+  std::complex<double> hifac = eifac*hffacta;
+
+  Hi[0] = hifac*std::sin(Theta)*std::sin(Phi);
+  Hi[1] = hifac*std::cos(Theta)*std::sin(Phi);
+  Hi[2] = hifac*std::cos(Phi);
+
+  for (i = 0; i < 3; i++) {
+    // electric field E [V m-1] = EF*E0
+    E[i] = Ei[i] + Es[i];
+    H[i] = Hi[i] + Hs[i];
+  }
+}
+
+// Calculate an - equation (5)
+std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
+                             std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                             std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
+
+  std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
+  std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
+
+  return Num/Denom;
+}
+
+// Calculate bn - equation (6)
+std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
+                             std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                             std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
+
+  std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
+  std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
+
+  return Num/Denom;
+}
+
+// Calculates S1 - equation (25a)
+std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
+                             double Pi, double Tau) {
+
+  return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
+}
+
+// Calculates S2 - equation (25b) (it's the same as (25a), just switches Pi and Tau)
+std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,
+                             double Pi, double Tau) {
+
+  return calc_S1(n, an, bn, Tau, Pi);
+}
+
+
+//**********************************************************************************//
+// This function calculates the Riccati-Bessel functions (Psi and Zeta) for a       //
+// real argument (x).                                                               //
+// Equations (20a) - (21b)                                                          //
+//                                                                                  //
+// Input parameters:                                                                //
+//   x: Real argument to evaluate Psi and Zeta                                      //
+//   nmax: Maximum number of terms to calculate Psi and Zeta                        //
+//                                                                                  //
+// Output parameters:                                                               //
+//   Psi, Zeta: Riccati-Bessel functions                                            //
+//**********************************************************************************//
+void calcPsiZeta(double x, int nmax,
+                 std::vector<std::complex<double> > D1,
+                 std::vector<std::complex<double> > D3,
+                 std::vector<std::complex<double> >& Psi,
+                 std::vector<std::complex<double> >& Zeta) {
+
+  int n;
+
+  //Upward recurrence for Psi and Zeta - equations (20a) - (21b)
+  Psi[0] = std::complex<double>(sin(x), 0);
+  Zeta[0] = std::complex<double>(sin(x), -cos(x));
+  for (n = 1; n <= nmax; n++) {
+    Psi[n] = Psi[n - 1]*(n/x - D1[n - 1]);
+    Zeta[n] = Zeta[n - 1]*(n/x - D3[n - 1]);
+  }
+}
+
+//**********************************************************************************//
+// This function calculates the logarithmic derivatives of the Riccati-Bessel       //
+// functions (D1 and D3) for a complex argument (z).                                //
+// Equations (16a), (16b) and (18a) - (18d)                                         //
+//                                                                                  //
+// Input parameters:                                                                //
+//   z: Complex argument to evaluate D1 and D3                                      //
+//   nmax: Maximum number of terms to calculate D1 and D3                           //
+//                                                                                  //
+// Output parameters:                                                               //
+//   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions                //
+//**********************************************************************************//
+void calcD1D3(std::complex<double> z, int nmax,
+              std::vector<std::complex<double> >& D1,
+              std::vector<std::complex<double> >& D3) {
+
+  int n;
+  std::complex<double> nz, PsiZeta;
+
+  // Downward recurrence for D1 - equations (16a) and (16b)
+  D1[nmax] = std::complex<double>(0.0, 0.0);
+  for (n = nmax; n > 0; n--) {
+    nz = double(n)/z;
+    D1[n - 1] = nz - 1.0/(D1[n] + nz);
+  }
+
+  // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
+  PsiZeta = 0.5*(1.0 - std::complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))*exp(-2.0*z.imag()));
+  D3[0] = std::complex<double>(0.0, 1.0);
+  for (n = 1; n <= nmax; n++) {
+    nz = double(n)/z;
+    PsiZeta = PsiZeta*(nz - D1[n - 1])*(nz - D3[n - 1]);
+    D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta;
+  }
+}
+
+//**********************************************************************************//
+// This function calculates Pi and Tau for all values of Theta.                     //
+// Equations (26a) - (26c)                                                          //
+//                                                                                  //
+// Input parameters:                                                                //
+//   nmax: Maximum number of terms to calculate Pi and Tau                          //
+//   nTheta: Number of scattering angles                                            //
+//   Theta: Array containing all the scattering angles where the scattering         //
+//          amplitudes will be calculated                                           //
+//                                                                                  //
+// Output parameters:                                                               //
+//   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
+//**********************************************************************************//
+void calcPiTau(int nmax, double Theta, std::vector<double>& Pi, std::vector<double>& Tau) {
+
+  int n;
+  //****************************************************//
+  // Equations (26a) - (26c)                            //
+  //****************************************************//
+  // Initialize Pi and Tau
+  Pi[0] = 1.0;
+  Tau[0] = cos(Theta);
+  // Calculate the actual values
+  if (nmax > 1) {
+    Pi[1] = 3*Tau[0]*Pi[0];
+    Tau[1] = 2*Tau[0]*Pi[1] - 3*Pi[0];
+    for (n = 2; n < nmax; n++) {
+      Pi[n] = ((n + n + 1)*Tau[0]*Pi[n - 1] - (n + 1)*Pi[n - 2])/n;
+      Tau[n] = (n + 1)*Tau[0]*Pi[n] - (n + 2)*Pi[n - 1];
+    }
+  }
+}
+
+//**********************************************************************************//
+// This function calculates 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(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
+                std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
+  //************************************************************************//
+  // Calculate the index of the first layer. It can be either 0 (default)   //
+  // or the index of the outermost PEC layer. In the latter case all layers //
+  // below the PEC are discarded.                                           //
+  //************************************************************************//
+
+  int fl = (pl > 0) ? pl : 0;
+
+  if (nmax <= 0) {
+    nmax = Nmax(L, fl, pl, x, m);
+  }
+
+  std::complex<double> z1, z2;
+  std::complex<double> Num, Denom;
+  std::complex<double> G1, G2;
+  std::complex<double> Temp;
+
+  int n, l;
+
+  //**************************************************************************//
+  // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which  //
+  // means that index = layer number - 1 or index = n - 1. The only exception //
+  // are the arrays for representing D1, D3 and Q because they need a value   //
+  // for the index 0 (zero), hence it is important to consider this shift     //
+  // between different arrays. The change was done to optimize memory usage.  //
+  //**************************************************************************//
+
+  // Allocate memory to the arrays
+  std::vector<std::vector<std::complex<double> > > D1_mlxl, D1_mlxlM1;
+  D1_mlxl.resize(L);
+  D1_mlxlM1.resize(L);
+
+  std::vector<std::vector<std::complex<double> > > D3_mlxl, D3_mlxlM1;
+  D3_mlxl.resize(L);
+  D3_mlxlM1.resize(L);
+
+  std::vector<std::vector<std::complex<double> > > Q;
+  Q.resize(L);
+
+  std::vector<std::vector<std::complex<double> > > Ha, Hb;
+  Ha.resize(L);
+  Hb.resize(L);
+
+  for (l = 0; l < L; l++) {
+    D1_mlxl[l].resize(nmax + 1);
+    D1_mlxlM1[l].resize(nmax + 1);
+
+    D3_mlxl[l].resize(nmax + 1);
+    D3_mlxlM1[l].resize(nmax + 1);
+
+    Q[l].resize(nmax + 1);
+
+    Ha[l].resize(nmax);
+    Hb[l].resize(nmax);
+  }
+
+  an.resize(nmax);
+  bn.resize(nmax);
+
+  std::vector<std::complex<double> > D1XL, D3XL;
+  D1XL.resize(nmax + 1);
+  D3XL.resize(nmax + 1);
+
+
+  std::vector<std::complex<double> > PsiXL, ZetaXL;
+  PsiXL.resize(nmax + 1);
+  ZetaXL.resize(nmax + 1);
+
+  //*************************************************//
+  // Calculate D1 and D3 for z1 in the first layer   //
+  //*************************************************//
+  if (fl == pl) {  // PEC layer
+    for (n = 0; n <= nmax; n++) {
+      D1_mlxl[fl][n] = std::complex<double>(0.0, -1.0);
+      D3_mlxl[fl][n] = std::complex<double>(0.0, 1.0);
+    }
+  } else { // Regular layer
+    z1 = x[fl]* m[fl];
+
+    // Calculate D1 and D3
+    calcD1D3(z1, nmax, D1_mlxl[fl], D3_mlxl[fl]);
+  }
+
+  //******************************************************************//
+  // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
+  //******************************************************************//
+  for (n = 0; n < nmax; n++) {
+    Ha[fl][n] = D1_mlxl[fl][n + 1];
+    Hb[fl][n] = D1_mlxl[fl][n + 1];
+  }
+
+  //*****************************************************//
+  // Iteration from the second layer to the last one (L) //
+  //*****************************************************//
+  for (l = fl + 1; l < L; l++) {
+    //************************************************************//
+    //Calculate D1 and D3 for z1 and z2 in the layers fl+1..L     //
+    //************************************************************//
+    z1 = x[l]*m[l];
+    z2 = x[l - 1]*m[l];
+
+    //Calculate D1 and D3 for z1
+    calcD1D3(z1, nmax, D1_mlxl[l], D3_mlxl[l]);
+
+    //Calculate D1 and D3 for z2
+    calcD1D3(z2, nmax, D1_mlxlM1[l], D3_mlxlM1[l]);
+
+    //*********************************************//
+    //Calculate Q, Ha and Hb in the layers fl+1..L //
+    //*********************************************//
+
+    // Upward recurrence for Q - equations (19a) and (19b)
+    Num = exp(-2.0*(z1.imag() - z2.imag()))*std::complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
+    Denom = std::complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
+    Q[l][0] = Num/Denom;
+
+    for (n = 1; n <= nmax; n++) {
+      Num = (z1*D1_mlxl[l][n] + double(n))*(double(n) - z1*D3_mlxl[l][n - 1]);
+      Denom = (z2*D1_mlxlM1[l][n] + double(n))*(double(n) - z2*D3_mlxlM1[l][n - 1]);
+
+      Q[l][n] = (((x[l - 1]*x[l - 1])/(x[l]*x[l])* Q[l][n - 1])*Num)/Denom;
+    }
+
+    // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
+    for (n = 1; n <= nmax; n++) {
+      //Ha
+      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+        G1 = -D1_mlxlM1[l][n];
+        G2 = -D3_mlxlM1[l][n];
+      } else {
+        G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[l][n]);
+        G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[l][n]);
+      }
+
+      Temp = Q[l][n]*G1;
+
+      Num = (G2*D1_mlxl[l][n]) - (Temp*D3_mlxl[l][n]);
+      Denom = G2 - Temp;
+
+      Ha[l][n - 1] = Num/Denom;
+
+      //Hb
+      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+        G1 = Hb[l - 1][n - 1];
+        G2 = Hb[l - 1][n - 1];
+      } else {
+        G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[l][n]);
+        G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[l][n]);
+      }
+
+      Temp = Q[l][n]*G1;
+
+      Num = (G2*D1_mlxl[l][n]) - (Temp* D3_mlxl[l][n]);
+      Denom = (G2- Temp);
+
+      Hb[l][n - 1] = (Num/ Denom);
+    }
+  }
+
+  //**************************************//
+  //Calculate D1, D3, Psi and Zeta for XL //
+  //**************************************//
+
+  // Calculate D1XL and D3XL
+  calcD1D3(x[L - 1], nmax, D1XL, D3XL);
+
+  // Calculate PsiXL and ZetaXL
+  calcPsiZeta(x[L - 1], nmax, D1XL, D3XL, PsiXL, ZetaXL);
+
+  //*********************************************************************//
+  // Finally, we calculate the scattering coefficients (an and bn) and   //
+  // the angular functions (Pi and Tau). Note that for these arrays the  //
+  // first layer is 0 (zero), in future versions all arrays will follow  //
+  // this convention to save memory. (13 Nov, 2014)                      //
+  //*********************************************************************//
+  for (n = 0; n < nmax; n++) {
+    //********************************************************************//
+    //Expressions for calculating an and bn coefficients are not valid if //
+    //there is only one PEC layer (ie, for a simple PEC sphere).          //
+    //********************************************************************//
+    if (pl < (L - 1)) {
+      an[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+      bn[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+    } else {
+      an[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+      bn[n] = PsiXL[n + 1]/ZetaXL[n + 1];
+    }
+  }
+
+  return nmax;
+}
+
+//**********************************************************************************//
+// This function calculates the actual scattering parameters and amplitudes         //
+//                                                                                  //
+// 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]     //
+//   nTheta: Number of scattering angles                                            //
+//   Theta: Array containing all the scattering angles where the scattering         //
+//          amplitudes will be calculated                                           //
+//   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:                                                               //
+//   Qext: Efficiency factor for extinction                                         //
+//   Qsca: Efficiency factor for scattering                                         //
+//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+//   Qbk: Efficiency factor for backscattering                                      //
+//   Qpr: Efficiency factor for the radiation pressure                              //
+//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+//   S1, S2: Complex scattering amplitudes                                          //
+//                                                                                  //
+// Return value:                                                                    //
+//   Number of multipolar expansion terms used for the calculations                 //
+//**********************************************************************************//
+
+int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
+         int nTheta, std::vector<double> Theta, 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 i, n, t;
+  std::vector<std::complex<double> > an, bn;
+  std::complex<double> Qbktmp;
+
+  // Calculate scattering coefficients
+  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
+
+  std::vector<double> Pi, Tau;
+  Pi.resize(nmax);
+  Tau.resize(nmax);
+
+  double x2 = x[L - 1]*x[L - 1];
+
+  // Initialize the scattering parameters
+  *Qext = 0;
+  *Qsca = 0;
+  *Qabs = 0;
+  *Qbk = 0;
+  Qbktmp = std::complex<double>(0.0, 0.0);
+  *Qpr = 0;
+  *g = 0;
+  *Albedo = 0;
+
+  // Initialize the scattering amplitudes
+  for (t = 0; t < nTheta; t++) {
+    S1[t] = std::complex<double>(0.0, 0.0);
+    S2[t] = std::complex<double>(0.0, 0.0);
+  }
+
+  // By using downward recurrence we avoid loss of precision due to float rounding errors
+  // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
+  //      http://en.wikipedia.org/wiki/Loss_of_significance
+  for (i = nmax - 2; i >= 0; i--) {
+    n = i + 1;
+    // Equation (27)
+    *Qext += (n + n + 1)*(an[i].real() + bn[i].real());
+    // Equation (28)
+    *Qsca += (n + n + 1)*(an[i].real()*an[i].real() + an[i].imag()*an[i].imag() + bn[i].real()*bn[i].real() + bn[i].imag()*bn[i].imag());
+    // Equation (29) TODO We must check carefully this equation. If we
+    // remove the typecast to double then the result changes. Which is
+    // the correct one??? Ovidio (2014/12/10) With cast ratio will
+    // give double, without cast (n + n + 1)/(n*(n + 1)) will be
+    // rounded to integer. Tig (2015/02/24)
+    *Qpr += ((n*(n + 2)/(n + 1))*((an[i]*std::conj(an[n]) + bn[i]*std::conj(bn[n])).real()) + ((double)(n + n + 1)/(n*(n + 1)))*(an[i]*std::conj(bn[i])).real());
+    // Equation (33)
+    Qbktmp = Qbktmp + (double)(n + n + 1)*(1 - 2*(n % 2))*(an[i]- bn[i]);
+
+    //****************************************************//
+    // Calculate the scattering amplitudes (S1 and S2)    //
+    // Equations (25a) - (25b)                            //
+    //****************************************************//
+    for (t = 0; t < nTheta; t++) {
+      calcPiTau(nmax, Theta[t], Pi, Tau);
+
+      S1[t] += calc_S1(n, an[i], bn[i], Pi[i], Tau[i]);
+      S2[t] += calc_S2(n, an[i], bn[i], Pi[i], Tau[i]);
+    }
+  }
+
+  *Qext = 2*(*Qext)/x2;                                 // Equation (27)
+  *Qsca = 2*(*Qsca)/x2;                                 // Equation (28)
+  *Qpr = *Qext - 4*(*Qpr)/x2;                           // Equation (29)
+
+  *Qabs = *Qext - *Qsca;                                // Equation (30)
+  *Albedo = *Qsca / *Qext;                              // Equation (31)
+  *g = (*Qext - *Qpr) / *Qsca;                          // Equation (32)
+
+  *Qbk = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2;    // Equation (33)
+
+  return nmax;
+}
+
+//**********************************************************************************//
+// This function is just a wrapper to call the full 'nMie' function with fewer      //
+// parameters, it is here mainly for compatibility with older versions of the       //
+// program. Also, you can use it if you neither have a PEC layer nor want to define //
+// any limit for the maximum number of terms.                                       //
+//                                                                                  //
+// Input parameters:                                                                //
+//   L: Number of layers                                                            //
+//   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]     //
+//   nTheta: Number of scattering angles                                            //
+//   Theta: Array containing all the scattering angles where the scattering         //
+//          amplitudes will be calculated                                           //
+//                                                                                  //
+// Output parameters:                                                               //
+//   Qext: Efficiency factor for extinction                                         //
+//   Qsca: Efficiency factor for scattering                                         //
+//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+//   Qbk: Efficiency factor for backscattering                                      //
+//   Qpr: Efficiency factor for the radiation pressure                              //
+//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+//   S1, S2: Complex scattering amplitudes                                          //
+//                                                                                  //
+// Return value:                                                                    //
+//   Number of multipolar expansion terms used for the calculations                 //
+//**********************************************************************************//
+
+int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
+         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) {
+
+  return nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+}
+
+
+//**********************************************************************************//
+// This function is just a wrapper to call the full 'nMie' function with fewer      //
+// parameters, it is useful if you want to include a PEC layer but not a limit      //
+// for the maximum number of terms.                                                 //
+//                                                                                  //
+// 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]     //
+//   nTheta: Number of scattering angles                                            //
+//   Theta: Array containing all the scattering angles where the scattering         //
+//          amplitudes will be calculated                                           //
+//                                                                                  //
+// Output parameters:                                                               //
+//   Qext: Efficiency factor for extinction                                         //
+//   Qsca: Efficiency factor for scattering                                         //
+//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+//   Qbk: Efficiency factor for backscattering                                      //
+//   Qpr: Efficiency factor for the radiation pressure                              //
+//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+//   S1, S2: Complex scattering amplitudes                                          //
+//                                                                                  //
+// Return value:                                                                    //
+//   Number of multipolar expansion terms used for the calculations                 //
+//**********************************************************************************//
+
+int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
+         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) {
+
+  return nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+}
+
+//**********************************************************************************//
+// This function is just a wrapper to call the full 'nMie' function with fewer      //
+// parameters, it is useful if you want to include a limit for the maximum number   //
+// of terms but not a PEC layer.                                                    //
+//                                                                                  //
+// Input parameters:                                                                //
+//   L: Number of layers                                                            //
+//   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]     //
+//   nTheta: Number of scattering angles                                            //
+//   Theta: Array containing all the scattering angles where the scattering         //
+//          amplitudes will be calculated                                           //
+//   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:                                                               //
+//   Qext: Efficiency factor for extinction                                         //
+//   Qsca: Efficiency factor for scattering                                         //
+//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+//   Qbk: Efficiency factor for backscattering                                      //
+//   Qpr: Efficiency factor for the radiation pressure                              //
+//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+//   S1, S2: Complex scattering amplitudes                                          //
+//                                                                                  //
+// Return value:                                                                    //
+//   Number of multipolar expansion terms used for the calculations                 //
+//**********************************************************************************//
+
+int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
+         int nTheta, std::vector<double> Theta, 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) {
+
+  return nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+}
+
+
+//**********************************************************************************//
+// This function calculates complex electric and magnetic field in the surroundings //
+// and inside (TODO) the particle.                                                  //
+//                                                                                  //
+// Input parameters:                                                                //
+//   L: Number of layers                                                            //
+//   pl: Index of PEC layer. If there is none just send 0 (zero)                    //
+//   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 0 (zero) and the function will calculate it.       //
+//   ncoord: Number of coordinate points                                            //
+//   Coords: Array containing all coordinates where the complex electric and        //
+//           magnetic fields will be calculated                                     //
+//                                                                                  //
+// Output parameters:                                                               //
+//   E, H: Complex electric and magnetic field at the provided coordinates          //
+//                                                                                  //
+// Return value:                                                                    //
+//   Number of multipolar expansion terms used for the calculations                 //
+//**********************************************************************************//
+
+int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax, int ncoord, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
+
+  int i, c;
+  double Rho, Phi, Theta;
+  std::vector<std::complex<double> > an, bn;
+
+  // This array contains the fields in spherical coordinates
+  std::vector<std::complex<double> > Es, Hs;
+  Es.resize(3);
+  Hs.resize(3);
+
+
+  // Calculate scattering coefficients
+  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
+
+  std::vector<double> Pi, Tau;
+  Pi.resize(nmax);
+  Tau.resize(nmax);
+
+  for (c = 0; c < ncoord; c++) {
+    // Convert to spherical coordinates
+    Rho = sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c] + Zp[c]*Zp[c]);
+
+    // Avoid convergence problems due to Rho too small
+    if (Rho < 1e-5) {
+      Rho = 1e-5;
+    }
+
+    //If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
+    if (Rho == 0.0) {
+      Theta = 0.0;
+    } else {
+      Theta = acos(Zp[c]/Rho);
+    }
+
+    //If Xp=Yp=0 then Phi is undefined. Just set it to zero to zero to avoid problems
+    if ((Xp[c] == 0.0) and (Yp[c] == 0.0)) {
+      Phi = 0.0;
+    } else {
+      Phi = acos(Xp[c]/sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c]));
+    }
+
+    calcPiTau(nmax, Theta, Pi, Tau);
+
+    //*******************************************************//
+    // external scattering field = incident + scattered      //
+    // BH p.92 (4.37), 94 (4.45), 95 (4.50)                  //
+    // assume: medium is non-absorbing; refim = 0; Uabs = 0  //
+    //*******************************************************//
+
+    // Firstly the easiest case: the field outside the particle
+    if (Rho >= x[L - 1]) {
+      fieldExt(nmax, Rho, Phi, Theta, Pi, Tau, an, bn, Es, Hs);
+    } else {
+      // TODO, for now just set all the fields to zero
+      for (i = 0; i < 3; i++) {
+        Es[i] = std::complex<double>(0.0, 0.0);
+        Hs[i] = std::complex<double>(0.0, 0.0);
+      }
+    }
+
+    //Now, convert the fields back to cartesian coordinates
+    E[c][0] = std::sin(Theta)*std::cos(Phi)*Es[0] + std::cos(Theta)*std::cos(Phi)*Es[1] - std::sin(Phi)*Es[2];
+    E[c][1] = std::sin(Theta)*std::sin(Phi)*Es[0] + std::cos(Theta)*std::sin(Phi)*Es[1] + std::cos(Phi)*Es[2];
+    E[c][2] = std::cos(Theta)*Es[0] - std::sin(Theta)*Es[1];
+
+    H[c][0] = std::sin(Theta)*std::cos(Phi)*Hs[0] + std::cos(Theta)*std::cos(Phi)*Hs[1] - std::sin(Phi)*Hs[2];
+    H[c][1] = std::sin(Theta)*std::sin(Phi)*Hs[0] + std::cos(Theta)*std::sin(Phi)*Hs[1] + std::cos(Phi)*Hs[2];
+    H[c][2] = std::cos(Theta)*Hs[0] - std::sin(Theta)*Hs[1];
+  }
+
+  return nmax;
+}

nmie-old.h

@@ -0,0 +1,58 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.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/>.         //
+//**********************************************************************************//
+
+#define VERSION "0.3.1"
+#include <complex>
+#include <vector>
+
+int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
+		        std::vector<std::complex<double> > &an, std::vector<std::complex<double> > &bn);
+
+int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
+         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(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
+         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(int L, std::vector<double> x, std::vector<std::complex<double> > m,
+         int nTheta, std::vector<double> Theta, 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(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
+         int nTheta, std::vector<double> Theta, 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 nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
+           int ncoord, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
+		   std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
+
+

nmie-wrapper.cc

@@ -71,6 +71,10 @@ namespace nmie {
       *Albedo = multi_layer_mie.GetAlbedo();
       S1 = multi_layer_mie.GetS1();
       S2 = multi_layer_mie.GetS2();
+      
+      printf("S1 = %16.14f + i*%16.14f, S1_ass =  %16.14f + i*%16.14f\n",
+             multi_layer_mie.GetS1()[0].real(), multi_layer_mie.GetS1()[0].imag(), S1[0].real(), S1[0].real());
+      
       //multi_layer_mie.GetFailed();
     } catch(const std::invalid_argument& ia) {
       // Will catch if  multi_layer_mie fails or other errors.
@@ -505,6 +509,7 @@ namespace nmie {
     }
     nmax_ += 15;  // Final nmax_ value
   }
+
   //**********************************************************************************//
   // This function calculates the spherical Bessel (jn) and Hankel (h1n) functions    //
   // and their derivatives for a given complex value z. See pag. 87 B&H.              //
@@ -796,7 +801,7 @@ c    MM + 1  and - 1, alternately
     MM = - 1; 
     KK = 2*N +3; //debug 3
 // c                                 ** Eq. R25b, k=2
-    CAK    = static_cast<std::complex<double> >(MM*KK) * ZINV; //debug -3 ZINV
+    CAK    = static_cast<std::complex<double> >(MM*KK)*ZINV; //debug -3 ZINV
     CDENOM = CAK;
     CNUMER = CDENOM + one/CONFRA; //-3zinv+z
     KOUNT  = 1;
@@ -807,15 +812,15 @@ c    MM + 1  and - 1, alternately
         throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
       }
       MM *= - 1;      KK += 2;  //debug  mm=1 kk=5
-      CAK = static_cast<std::complex<double> >(MM*KK) * ZINV; //    ** Eq. R25b //debug 5zinv
+      CAK = static_cast<std::complex<double> >(MM*KK)*ZINV; //    ** Eq. R25b //debug 5zinv
      //  //c ** Eq. R32    Ill-conditioned case -- stride two terms instead of one
      //  if (std::abs(CNUMER/CAK) >= EPS1 ||  std::abs(CDENOM/CAK) >= EPS1) {
      //         //c                       ** Eq. R34
-     //         CNTN   = CAK * CNUMER + 1.0;
-     //         CDTD   = CAK * CDENOM + 1.0;
-     //         CONFRA = (CNTN/CDTD) * CONFRA; // ** Eq. R33
+     //         CNTN   = CAK*CNUMER + 1.0;
+     //         CDTD   = CAK*CDENOM + 1.0;
+     //         CONFRA = (CNTN/CDTD)*CONFRA; // ** Eq. R33
      //         MM  *= - 1;        KK  += 2;
-     //         CAK = static_cast<std::complex<double> >(MM*KK) * ZINV; // ** Eq. R25b
+     //         CAK = static_cast<std::complex<double> >(MM*KK)*ZINV; // ** Eq. R25b
      //         //c                        ** Eq. R35
      //         CNUMER = CAK + CNUMER/CNTN;
      //         CDENOM = CAK + CDENOM/CDTD;
@@ -826,7 +831,7 @@ c    MM + 1  and - 1, alternately
       {
         CAPT   = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
         // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
-       CONFRA = CAPT * CONFRA; // ** Eq. R26
+       CONFRA = CAPT*CONFRA; // ** Eq. R26
        //if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
        //c                                  ** Check for convergence; Eq. R31
        if (std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2) {
@@ -931,6 +936,7 @@ c    MM + 1  and - 1, alternately
       //calcSinglePiTau(std::cos(theta_[t]), Pi[t], Tau[t]); // It is slow!!
     }
   }  // end of void MultiLayerMie::calcAllPiTau(...)
+
   //**********************************************************************************//
   // This function calculates the scattering coefficients required to calculate       //
   // both the near- and far-field parameters.                                         //
@@ -950,7 +956,7 @@ c    MM + 1  and - 1, alternately
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  void MultiLayerMie::ScattCoeffs(std::vector<std::complex<double> >& an,
+  void MultiLayerMie::ExtScattCoeffs(std::vector<std::complex<double> >& an,
                                   std::vector<std::complex<double> >& bn) {
     const std::vector<double>& x = size_parameter_;
     const std::vector<std::complex<double> >& m = index_;
@@ -966,6 +972,8 @@ c    MM + 1  and - 1, alternately
     // int fl = (pl > - 1) ? pl : 0;
     // This will give the same result, however, it corresponds the
     // logic - if there is PEC, than first layer is PEC.
+    // Well, I followed the logic: First layer is always zero unless it has 
+    // an upper PEC layer.
     int fl = (pl > 0) ? pl : 0;
     if (nmax_ <= 0) Nmax(fl);
 
@@ -979,22 +987,23 @@ c    MM + 1  and - 1, alternately
     //**************************************************************************//
     // Allocate memory to the arrays
     std::vector<std::complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
-      D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
+                                       D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
+
     std::vector<std::vector<std::complex<double> > > Q(L), Ha(L), Hb(L);
+
     for (int l = 0; l < L; l++) {
-      // D1_mlxl[l].resize(nmax_ + 1);
-      // D1_mlxlM1[l].resize(nmax_ + 1);
-      // D3_mlxl[l].resize(nmax_ + 1);
-      // D3_mlxlM1[l].resize(nmax_ + 1);
       Q[l].resize(nmax_ + 1);
       Ha[l].resize(nmax_);
       Hb[l].resize(nmax_);
     }
+
     an.resize(nmax_);
     bn.resize(nmax_);
     PsiZeta_.resize(nmax_ + 1);
+
     std::vector<std::complex<double> > D1XL(nmax_ + 1), D3XL(nmax_ + 1), 
-      PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
+                                       PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
+
     //*************************************************//
     // Calculate D1 and D3 for z1 in the first layer   //
     //*************************************************//
@@ -1044,7 +1053,7 @@ c    MM + 1  and - 1, alternately
       //*********************************************//
       // Upward recurrence for Q - equations (19a) and (19b)
       Num = std::exp(-2.0*(z1.imag() - z2.imag()))
-        * std::complex<double>(std::cos(-2.0*z2.real()) - std::exp(-2.0*z2.imag()), std::sin(-2.0*z2.real()));
+       *std::complex<double>(std::cos(-2.0*z2.real()) - std::exp(-2.0*z2.imag()), std::sin(-2.0*z2.real()));
       Denom = std::complex<double>(std::cos(-2.0*z1.real()) - std::exp(-2.0*z1.imag()), std::sin(-2.0*z1.real()));
       Q[l][0] = Num/Denom;
       for (int n = 1; n <= nmax_; n++) {
@@ -1108,7 +1117,7 @@ c    MM + 1  and - 1, alternately
         bn[n] = PsiXL[n + 1]/ZetaXL[n + 1];
       }
     }  // end of for an and bn terms
-  }  // end of void MultiLayerMie::ScattCoeffs(...)
+  }  // end of void MultiLayerMie::ExtScattCoeffs(...)
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
@@ -1138,7 +1147,7 @@ c    MM + 1  and - 1, alternately
     Qbk_ch_norm_.resize(nmax_ - 1);
     Qpr_ch_norm_.resize(nmax_ - 1);
     // Initialize the scattering amplitudes
-    std::vector<std::complex<double> >        tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
+    std::vector<std::complex<double> > tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
     S1_.swap(tmp1);
     S2_ = S1_;
   }
@@ -1162,15 +1171,15 @@ c    MM + 1  and - 1, alternately
   //                                                                                  //
   // 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]     //
+  //   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]     //
   //   nTheta: Number of scattering angles                                            //
   //   Theta: Array containing all the scattering angles where the scattering         //
   //          amplitudes will be calculated                                           //
-  //   nmax_: Maximum number of multipolar expansion terms to be used for the          //
+  //   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              //
+  //         set this parameter to -1 and the function will calculate it              //
   //                                                                                  //
   // Output parameters:                                                               //
   //   Qext: Efficiency factor for extinction                                         //
@@ -1195,7 +1204,7 @@ c    MM + 1  and - 1, alternately
       throw std::invalid_argument("Initialize model first!");
     const std::vector<double>& x = size_parameter_;
     // Calculate scattering coefficients
-    ScattCoeffs(an_, bn_);
+    ExtScattCoeffs(an_, bn_);
 
     // std::vector< std::vector<double> > Pi(nmax_), Tau(nmax_);
     std::vector< std::vector<double> > Pi, Tau;
@@ -1276,7 +1285,7 @@ c    MM + 1  and - 1, alternately
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void MultiLayerMie::ScattCoeffsLayerdInit() {
+  void MultiLayerMie::IntScattCoeffsInit() {
     const int L = index_.size();
     // we need to fill
     // std::vector< std::vector<std::complex<double> > > al_n_, bl_n_, cl_n_, dl_n_;
@@ -1300,25 +1309,25 @@ c    MM + 1  and - 1, alternately
       bl_n_[L][i] = bn_[i];
       cl_n_[L][i] = c_one;
       dl_n_[L][i] = c_one;
-      if (i<3) printf(" (%g) ", std::abs(an_[i]));
+      if (i < 3) printf(" (%g) ", std::abs(an_[i]));
     }
 
   }
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void MultiLayerMie::ScattCoeffsLayerd() {
+  void MultiLayerMie::IntScattCoeffs() {
     if (!isMieCalculated_)
-      throw std::invalid_argument("(ScattCoeffsLayerd) You should run calculations first!");
-    ScattCoeffsLayerdInit();
+      throw std::invalid_argument("(IntScattCoeffs) You should run calculations first!");
+    IntScattCoeffsInit();
     const int L = index_.size();
     std::vector<std::complex<double> > z(L), z1(L);
     for (int i = 0; i < L - 1; ++i) {
       z[i]  =size_parameter_[i]*index_[i];
       z1[i]=size_parameter_[i]*index_[i + 1];
     }
-    z[L - 1]  =size_parameter_[L - 1]*index_[L - 1];
-    z1[L - 1]  =size_parameter_[L - 1];
+    z[L - 1] = size_parameter_[L - 1]*index_[L - 1];
+    z1[L - 1] = size_parameter_[L - 1];
     std::vector< std::vector<std::complex<double> > > D1z(L), D1z1(L), D3z(L), D3z1(L);
     std::vector< std::vector<std::complex<double> > > Psiz(L), Psiz1(L), Zetaz(L), Zetaz1(L);
     for (int l = 0; l < L; ++l) {
@@ -1345,37 +1354,31 @@ c    MM + 1  and - 1, alternately
     for (int l = L - 1; l >= 0; --l) {
       for (int n = 0; n < nmax_; ++n) {
         // al_n
-        auto denom = m1[l]*Zetaz[l][n + 1] * (D1z[l][n + 1] - D3z[l][n + 1]);
-        al_n_[l][n] = D1z[l][n + 1]* m1[l]
-          *(al_n_[l + 1][n]*Zetaz1[l][n + 1] - dl_n_[l + 1][n]*Psiz1[l][n + 1])
-          - m[l]*(-D1z1[l][n + 1]*dl_n_[l + 1][n]*Psiz1[l][n + 1]
-                  +D3z1[l][n + 1]*al_n_[l + 1][n]*Zetaz1[l][n + 1]);
+        auto denom = m1[l]*Zetaz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+        al_n_[l][n] = D1z[l][n + 1]*m1[l]*(al_n_[l + 1][n]*Zetaz1[l][n + 1] - dl_n_[l + 1][n]*Psiz1[l][n + 1])
+                      - m[l]*(-D1z1[l][n + 1]*dl_n_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*al_n_[l + 1][n]*Zetaz1[l][n + 1]);
         al_n_[l][n] /= denom;
-        // if (n<2) printf("denom[%d][%d]:%g \n", l, n,
-        //                   std::abs(Psiz[l][n + 1]));
+
         // dl_n
-        denom = m1[l]*Psiz[l][n + 1] * (D1z[l][n + 1] - D3z[l][n + 1]);
-        dl_n_[l][n] = D3z[l][n + 1]*m1[l]
-          *(al_n_[l + 1][n]*Zetaz1[l][n + 1] - dl_n_[l + 1][n]*Psiz1[l][n + 1])
-          - m[l]*(-D1z1[l][n + 1]*dl_n_[l + 1][n]*Psiz1[l][n + 1]
-                  +D3z1[l][n + 1]*al_n_[l + 1][n]*Zetaz1[l][n + 1]);
+        denom = m1[l]*Psiz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+        dl_n_[l][n] = D3z[l][n + 1]*m1[l]*(al_n_[l + 1][n]*Zetaz1[l][n + 1] - dl_n_[l + 1][n]*Psiz1[l][n + 1])
+                      - m[l]*(-D1z1[l][n + 1]*dl_n_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*al_n_[l + 1][n]*Zetaz1[l][n + 1]);
         dl_n_[l][n] /= denom;
+
         // bl_n
-        denom = m1[l]*Zetaz[l][n + 1] * (D1z[l][n + 1] - D3z[l][n + 1]);
-        bl_n_[l][n] = D1z[l][n + 1]* m[l]
-          *(bl_n_[l + 1][n]*Zetaz1[l][n + 1] - cl_n_[l + 1][n]*Psiz1[l][n + 1])
-          - m1[l]*(-D1z1[l][n + 1]*cl_n_[l + 1][n]*Psiz1[l][n + 1]
-                  +D3z1[l][n + 1]*bl_n_[l + 1][n]*Zetaz1[l][n + 1]);
+        denom = m1[l]*Zetaz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+        bl_n_[l][n] = D1z[l][n + 1]*m[l]*(bl_n_[l + 1][n]*Zetaz1[l][n + 1] - cl_n_[l + 1][n]*Psiz1[l][n + 1])
+                      - m1[l]*(-D1z1[l][n + 1]*cl_n_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*bl_n_[l + 1][n]*Zetaz1[l][n + 1]);
         bl_n_[l][n] /= denom;
+
         // cl_n
-        denom = m1[l]*Psiz[l][n + 1] * (D1z[l][n + 1] - D3z[l][n + 1]);
-        cl_n_[l][n] = D3z[l][n + 1]*m[l]
-          *(bl_n_[l + 1][n]*Zetaz1[l][n + 1] - cl_n_[l + 1][n]*Psiz1[l][n + 1])
-          - m1[l]*(-D1z1[l][n + 1]*cl_n_[l + 1][n]*Psiz1[l][n + 1]
-                  +D3z1[l][n + 1]*bl_n_[l + 1][n]*Zetaz1[l][n + 1]);
+        denom = m1[l]*Psiz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+        cl_n_[l][n] = D3z[l][n + 1]*m[l]*(bl_n_[l + 1][n]*Zetaz1[l][n + 1] - cl_n_[l + 1][n]*Psiz1[l][n + 1])
+                      - m1[l]*(-D1z1[l][n + 1]*cl_n_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*bl_n_[l + 1][n]*Zetaz1[l][n + 1]);
         cl_n_[l][n] /= denom;   
       }  // end of all n
     }  // end of for all l
+
     // Check the result and change  an__0 and bn__0 for exact zero
     for (int n = 0; n < nmax_; ++n) {
       if (std::abs(al_n_[0][n]) < 1e-10) al_n_[0][n] = 0.0;
@@ -1597,10 +1600,8 @@ c    MM + 1  and - 1, alternately
       for (int i = 0; i < 3; i++) {
         // if (n<3 && i==0) printf("\nn=%d",n);
         // if (n<3) printf("\nbefore !El[%d]=%g,%g! ", i, El[i].real(), El[i].imag());
-        Ei[i] = encap*(
-                       cl_n_[l][n]*vm1o1n[i] - c_i*dl_n_[l][n]*vn1e1n[i]
-          + c_i*al_n_[l][n]*vn3e1n[i] - bl_n_[l][n]*vm3o1n[i]
-);
+        Ei[i] = encap*(cl_n_[l][n]*vm1o1n[i] - c_i*dl_n_[l][n]*vn1e1n[i]
+                       + c_i*al_n_[l][n]*vn3e1n[i] - bl_n_[l][n]*vm3o1n[i]);
         El[i] = El[i] + encap*(cl_n_[l][n]*vm1o1n[i] - c_i*dl_n_[l][n]*vn1e1n[i]
                                + c_i*al_n_[l][n]*vn3e1n[i] - bl_n_[l][n]*vm3o1n[i]);
         Hl[i] = Hl[i] + encap*(-dl_n_[l][n]*vm1e1n[i] - c_i*cl_n_[l][n]*vn1o1n[i]
@@ -1664,7 +1665,7 @@ c    MM + 1  and - 1, alternately
     // Calculate scattering coefficients an_ and bn_
     RunMieCalculations();
     //nmax_=10;
-    ScattCoeffsLayerd();
+    IntScattCoeffs();
 
     std::vector<double> Pi(nmax_), Tau(nmax_);
     long total_points = coords_sp_[0].size();

nmie-wrapper.h

@@ -120,8 +120,7 @@ namespace nmie {
       GetFieldE(){return E_field_;};   // {X[], Y[], Z[]}
     std::vector<std::vector< std::complex<double> > >
       GetFieldH(){return H_field_;};
-    std::vector< std::vector<double> >   GetSpectra(double from_WL, double to_WL,
-                                                   int samples);  // ext, sca, abs, bk
+    std::vector< std::vector<double> > GetSpectra(double from_WL, double to_WL, int samples);  // ext, sca, abs, bk
     double GetRCSext();
     double GetRCSsca();
     double GetRCSabs();
@@ -129,19 +128,16 @@ namespace nmie {
     std::vector<double> GetPatternEk();
     std::vector<double> GetPatternHk();
     std::vector<double> GetPatternUnpolarized();
-    
-
 
     // Size parameter units
-    std::vector<double>                  GetLayerWidthSP();
+    std::vector<double> GetLayerWidthSP();
     // Same as to get target and coating index
-    std::vector< std::complex<double> >  GetLayerIndex();  
-    std::vector< std::array<double,3> >   GetFieldPointsSP();
+    std::vector< std::complex<double> > GetLayerIndex();  
+    std::vector< std::array<double,3> > GetFieldPointsSP();
     // Do we need normalize field to size parameter?
     /* std::vector<std::vector<std::complex<double> > >  GetFieldESP(); */
     /* std::vector<std::vector<std::complex<double> > >  GetFieldHSP(); */
-    std::vector< std::array<double,5> >   GetSpectraSP(double from_SP, double to_SP,
-						       int samples);  // WL,ext, sca, abs, bk
+    std::vector< std::array<double,5> > GetSpectraSP(double from_SP, double to_SP, int samples);  // WL,ext, sca, abs, bk
     double GetQext();
     double GetQsca();
     double GetQabs();
@@ -210,14 +206,16 @@ namespace nmie {
 			             std::vector<double>& Tau);
     void calcAllPiTau(std::vector< std::vector<double> >& Pi,
 		              std::vector< std::vector<double> >& Tau);
-    void ScattCoeffs(std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn); 
-    void ScattCoeffsLayerd();
-    void ScattCoeffsLayerdInit();
+    void ExtScattCoeffs(std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn); 
+    void IntScattCoeffs();
+    void IntScattCoeffsInit();
 
     void fieldExt(const double Rho, const double Phi, const double Theta, const  std::vector<double>& Pi, const std::vector<double>& Tau, std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
 
     void fieldInt(const double Rho, const double Phi, const double Theta, const  std::vector<double>& Pi, const std::vector<double>& Tau, std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
     
+    bool areIntCoeffsCalc_ = false;
+    bool areExtCoeffsCalc_ = false;
     bool isMieCalculated_ = false;
     double wavelength_ = 1.0;
     double total_radius_ = 0.0;
@@ -232,7 +230,7 @@ namespace nmie {
     std::vector<double> theta_;
     // Should be -1 if there is no PEC.
     int PEC_layer_position_ = -1;
-    // Set nmax_ manualy with SetMaxTermsNumber(int nmax) or in ScattCoeffs(..)
+    // Set nmax_ manualy with SetMaxTermsNumber(int nmax) or in ExtScattCoeffs(..)
     // with Nmax(int first_layer);
     int nmax_ = -1;
     int nmax_used_ = -1;

nmie.cc

@@ -1,5 +1,6 @@
 //**********************************************************************************//
 //    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
 //                                                                                  //
 //    This file is part of scattnlay                                                //
 //                                                                                  //
@@ -25,9 +26,9 @@
 //**********************************************************************************//
 
 //**********************************************************************************//
-// This library implements the algorithm for a multilayered sphere described by:    //
+// This class implements the algorithm for a multilayered sphere described by:      //
 //    [1] W. Yang, "Improved recursive algorithm for light scattering by a          //
-//        multilayered sphere,” Applied Optics,  vol. 42, Mar. 2003, pp. 1710-1720. //
+//        multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp. 1710-1720.  //
 //                                                                                  //
 // You can find the description of all the used equations in:                       //
 //    [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
@@ -36,970 +37,1500 @@
 //                                                                                  //
 // Hereinafter all equations numbers refer to [2]                                   //
 //**********************************************************************************//
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
 #include "nmie.h"
+#include <array>
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <stdexcept>
+#include <vector>
+
+namespace nmie {
+  //helpers
+  template<class T> inline T pow2(const T value) {return value*value;}
+
+  int round(double x) {
+    return x >= 0 ? (int)(x + 0.5):(int)(x - 0.5);
+  }
 
-#define round(x) ((x) >= 0 ? (int)((x) + 0.5):(int)((x) - 0.5))
 
-const double PI=3.14159265358979323846;
-// light speed [m s-1]
-double const cc = 2.99792458e8;
-// assume non-magnetic (MU=MU0=const) [N A-2]
-double const mu = 4.0*PI*1.0e-7;
+  //**********************************************************************************//
+  // This function emulates a C call to calculate the actual scattering parameters    //
+  // and amplitudes.                                                                  //
+  //                                                                                  //
+  // 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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //   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:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  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, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+
+    if (x.size() != L || m.size() != L)
+        throw std::invalid_argument("Declared number of layers do not fit x and m!");
+    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.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();
+    } catch(const std::invalid_argument& ia) {
+      // Will catch if  multi_layer_mie fails or other errors.
+      std::cerr << "Invalid argument: " << ia.what() << std::endl;
+      throw std::invalid_argument(ia);
+      return -1;
+    }
 
-// Calculate Nstop - equation (17)
-int Nstop(double xL) {
-  int result;
+    return 0;
+  }
 
-  if (xL <= 8) {
-    result = round(xL + 4*pow(xL, 1.0/3.0) + 1);
-  } else if (xL <= 4200) {
-    result = round(xL + 4.05*pow(xL, 1.0/3.0) + 2);
-  } else {
-    result = round(xL + 4*pow(xL, 1.0/3.0) + 2);
+  //**********************************************************************************//
+  // This function is just a wrapper to call the full 'nMie' function with fewer      //
+  // parameters, it is here mainly for compatibility with older versions of the       //
+  // program. Also, you can use it if you neither have a PEC layer nor want to define //
+  // any limit for the maximum number of terms.                                       //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  int nMie(const int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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) {
+    return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
-  return result;
-}
 
-//**********************************************************************************//
-int Nmax(int L, int fl, int pl,
-         std::vector<double> x,
-         std::vector<std::complex<double> > m) {
-  int i, result, ri, riM1;
-  result = Nstop(x[L - 1]);
-  for (i = fl; i < L; i++) {
-    if (i > pl) {
-      ri = round(std::abs(x[i]*m[i]));
-    } else {
-      ri = 0;
-    }
-    if (result < ri) {
-      result = ri;
-    }
+  //**********************************************************************************//
+  // This function is just a wrapper to call the full 'nMie' function with fewer      //
+  // parameters, it is useful if you want to include a PEC layer but not a limit      //
+  // for the maximum number of terms.                                                 //
+  //                                                                                  //
+  // 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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  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, 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) {
+    return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+  }
 
-    if ((i > fl) && ((i - 1) > pl)) {
-      riM1 = round(std::abs(x[i - 1]* m[i]));
-    } else {
-      riM1 = 0;
-    }
-    if (result < riM1) {
-      result = riM1;
+  //**********************************************************************************//
+  // This function is just a wrapper to call the full 'nMie' function with fewer      //
+  // parameters, it is useful if you want to include a limit for the maximum number   //
+  // of terms but not a PEC layer.                                                    //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //   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:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  int nMie(const int L, 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, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+    return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+  }
+
+
+  //**********************************************************************************//
+  // This function emulates a C call to calculate complex electric and magnetic field //
+  // in the surroundings and inside (TODO) the particle.                              //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   pl: Index of PEC layer. If there is none just send 0 (zero)                    //
+  //   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 0 (zero) and the function will calculate it.       //
+  //   ncoord: Number of coordinate points                                            //
+  //   Coords: Array containing all coordinates where the complex electric and        //
+  //           magnetic fields will be calculated                                     //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   E, H: Complex electric and magnetic field at the provided coordinates          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  int nField(const int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const int ncoord, const std::vector<double>& Xp_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
+    if (x.size() != L || m.size() != L)
+      throw std::invalid_argument("Declared number of layers do not fit x and m!");
+    if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
+        || E.size() != ncoord || H.size() != ncoord)
+      throw std::invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
+    for (auto f:E)
+      if (f.size() != 3)
+        throw std::invalid_argument("Field E is not 3D!");
+    for (auto f:H)
+      if (f.size() != 3)
+        throw std::invalid_argument("Field H is not 3D!");
+    try {
+      MultiLayerMie multi_layer_mie;
+      //multi_layer_mie.SetPECLayer(pl);
+      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();
+      //multi_layer_mie.GetFailed();
+    } catch(const std::invalid_argument& ia) {
+      // Will catch if  multi_layer_mie fails or other errors.
+      std::cerr << "Invalid argument: " << ia.what() << std::endl;
+      throw std::invalid_argument(ia);
+      return - 1;
     }
+
+    return 0;
   }
-  return result + 15;
-}
 
-//**********************************************************************************//
-// This function calculates the spherical Bessel (jn) and Hankel (h1n) functions    //
-// and their derivatives for a given complex value z. See pag. 87 B&H.              //
-//                                                                                  //
-// Input parameters:                                                                //
-//   z: Real argument to evaluate jn and h1n                                        //
-//   nmax: Maximum number of terms to calculate jn and h1n                          //
-//                                                                                  //
-// Output parameters:                                                               //
-//   jn, h1n: Spherical Bessel and Hankel functions                                 //
-//   jnp, h1np: Derivatives of the spherical Bessel and Hankel functions            //
-//                                                                                  //
-// The implementation follows the algorithm by I.J. Thompson and A.R. Barnett,      //
-// Comp. Phys. Comm. 47 (1987) 245-257.                                             //
-//                                                                                  //
-// Complex spherical Bessel functions from n=0..nmax-1 for z in the upper half      //
-// plane (Im(z) > -3).                                                              //
-//                                                                                  //
-//     j[n]   = j/n(z)                Regular solution: j[0]=sin(z)/z               //
-//     j'[n]  = d[j/n(z)]/dz                                                        //
-//     h1[n]  = h[0]/n(z)             Irregular Hankel function:                    //
-//     h1'[n] = d[h[0]/n(z)]/dz                h1[0] = j0(z) + i*y0(z)              //
-//                                                   = (sin(z)-i*cos(z))/z          //
-//                                                   = -i*exp(i*z)/z                //
-// Using complex CF1, and trigonometric forms for n=0 solutions.                    //
-//**********************************************************************************//
-int sbesjh(std::complex<double> z, int nmax, std::vector<std::complex<double> >& jn, std::vector<std::complex<double> >& jnp, std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np) {
 
-  const int limit = 20000;
-  double const accur = 1.0e-12;
-  double const tm30 = 1e-30;
+  // ********************************************************************** //
+  // Returns previously calculated Qext                                     //
+  // ********************************************************************** //
+  double MultiLayerMie::GetQext() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return Qext_;
+  }
 
-  int n;
-  double absc;
-  std::complex<double> zi, w;
-  std::complex<double> pl, f, b, d, c, del, jn0, jndb, h1nldb, h1nbdb;
 
-  absc = std::abs(std::real(z)) + std::abs(std::imag(z));
-  if ((absc < accur) || (std::imag(z) < -3.0)) {
-    return -1;
+  // ********************************************************************** //
+  // Returns previously calculated Qabs                                     //
+  // ********************************************************************** //
+  double MultiLayerMie::GetQabs() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return Qabs_;
   }
 
-  zi = 1.0/z;
-  w = zi + zi;
 
-  pl = double(nmax)*zi;
+  // ********************************************************************** //
+  // Returns previously calculated Qsca                                     //
+  // ********************************************************************** //
+  double MultiLayerMie::GetQsca() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return Qsca_;
+  }
 
-  f = pl + zi;
-  b = f + f + zi;
-  d = 0.0;
-  c = f;
-  for (n = 0; n < limit; n++) {
-    d = b - d;
-    c = b - 1.0/c;
 
-    absc = std::abs(std::real(d)) + std::abs(std::imag(d));
-    if (absc < tm30) {
-      d = tm30;
-    }
+  // ********************************************************************** //
+  // Returns previously calculated Qbk                                      //
+  // ********************************************************************** //
+  double MultiLayerMie::GetQbk() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return Qbk_;
+  }
 
-    absc = std::abs(std::real(c)) + std::abs(std::imag(c));
-    if (absc < tm30) {
-      c = tm30;
-    }
 
-    d = 1.0/d;
-    del = d*c;
-    f = f*del;
-    b += w;
+  // ********************************************************************** //
+  // Returns previously calculated Qpr                                      //
+  // ********************************************************************** //
+  double MultiLayerMie::GetQpr() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return Qpr_;
+  }
 
-    absc = std::abs(std::real(del - 1.0)) + std::abs(std::imag(del - 1.0));
 
-    if (absc < accur) {
-      // We have obtained the desired accuracy
-      break;
-    }
+  // ********************************************************************** //
+  // Returns previously calculated assymetry factor                         //
+  // ********************************************************************** //
+  double MultiLayerMie::GetAsymmetryFactor() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return asymmetry_factor_;
   }
 
-  if (absc > accur) {
-    // We were not able to obtain the desired accuracy
-    return -2;
+
+  // ********************************************************************** //
+  // Returns previously calculated Albedo                                   //
+  // ********************************************************************** //
+  double MultiLayerMie::GetAlbedo() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return albedo_;
   }
 
-  jn[nmax - 1] = tm30;
-  jnp[nmax - 1] = f*jn[nmax - 1];
 
-  // Downward recursion to n=0 (N.B.  Coulomb Functions)
-  for (n = nmax - 2; n >= 0; n--) {
-    jn[n] = pl*jn[n + 1] + jnp[n + 1];
-    jnp[n] = pl*jn[n] - jn[n + 1];
-    pl = pl - zi;
+  // ********************************************************************** //
+  // Returns previously calculated S1                                       //
+  // ********************************************************************** //
+  std::vector<std::complex<double> > MultiLayerMie::GetS1() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return S1_;
   }
 
-  // Calculate the n=0 Bessel Functions
-  jn0 = zi*std::sin(z);
-  h1n[0] = std::exp(std::complex<double>(0.0, 1.0)*z)*zi*(-std::complex<double>(0.0, 1.0));
-  h1np[0] = h1n[0]*(std::complex<double>(0.0, 1.0) - zi);
 
-  // Rescale j[n], j'[n], converting to spherical Bessel functions.
-  // Recur   h1[n], h1'[n] as spherical Bessel functions.
-  w = 1.0/jn[0];
-  pl = zi;
-  for (n = 0; n < nmax; n++) {
-    jn[n] = jn0*(w*jn[n]);
-    jnp[n] = jn0*(w*jnp[n]) - zi*jn[n];
-    if (n != 0) {
-      h1n[n] = (pl - zi)*h1n[n - 1] - h1np[n - 1];
+  // ********************************************************************** //
+  // Returns previously calculated S2                                       //
+  // ********************************************************************** //
+  std::vector<std::complex<double> > MultiLayerMie::GetS2() {
+    if (!isMieCalculated_)
+      throw std::invalid_argument("You should run calculations before result request!");
+    return S2_;
+  }
 
-      // check if hankel is increasing (upward stable)
-      if (std::abs(h1n[n]) < std::abs(h1n[n - 1])) {
-        jndb = z;
-        h1nldb = h1n[n];
-        h1nbdb = h1n[n - 1];
-      }
 
-      pl += zi;
+  // ********************************************************************** //
+  // Modify scattering (theta) angles                                       //
+  // ********************************************************************** //
+  void MultiLayerMie::SetAngles(const std::vector<double>& angles) {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    theta_ = angles;
+  }
 
-      h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
+
+  // ********************************************************************** //
+  // Modify size of all layers                                             //
+  // ********************************************************************** //
+  void MultiLayerMie::SetLayersSize(const std::vector<double>& layer_size) {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    size_param_.clear();
+    double prev_layer_size = 0.0;
+    for (auto curr_layer_size : layer_size) {
+      if (curr_layer_size <= 0.0)
+        throw std::invalid_argument("Size parameter should be positive!");
+      if (prev_layer_size > curr_layer_size)
+        throw std::invalid_argument
+          ("Size parameter for next layer should be larger than the previous one!");
+      prev_layer_size = curr_layer_size;
+      size_param_.push_back(curr_layer_size);
     }
   }
 
-  // success
-  return 0;
-}
-
-//**********************************************************************************//
-// This function calculates the spherical Bessel functions (bj and by) and the      //
-// logarithmic derivative (bd) for a given complex value z. See pag. 87 B&H.        //
-//                                                                                  //
-// Input parameters:                                                                //
-//   z: Complex argument to evaluate bj, by and bd                                  //
-//   nmax: Maximum number of terms to calculate bj, by and bd                       //
-//                                                                                  //
-// Output parameters:                                                               //
-//   bj, by: Spherical Bessel functions                                             //
-//   bd: Logarithmic derivative                                                     //
-//**********************************************************************************//
-void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<double> >& bj, std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd) {
-
-    std::vector<std::complex<double> > jn, jnp, h1n, h1np;
-    jn.resize(nmax);
-    jnp.resize(nmax);
-    h1n.resize(nmax);
-    h1np.resize(nmax);
-
-    // TODO verify that the function succeeds
-    int ifail = sbesjh(z, nmax, jn, jnp, h1n, h1np);
-
-    for (int n = 0; n < nmax; n++) {
-      bj[n] = jn[n];
-      by[n] = (h1n[n] - jn[n])/std::complex<double>(0.0, 1.0);
-      bd[n] = jnp[n]/jn[n] + 1.0/z;
-    }
-}
-
-// external scattering field = incident + scattered
-// BH p.92 (4.37), 94 (4.45), 95 (4.50)
-// assume: medium is non-absorbing; refim = 0; Uabs = 0
-void fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double> Pi, std::vector<double> Tau,
-              std::vector<std::complex<double> > an, std::vector<std::complex<double> > bn,
-              std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
-
-  int i, n, n1;
-  double rn;
-  std::complex<double> ci, zn, xxip, encap;
-  std::vector<std::complex<double> > vm3o1n, vm3e1n, vn3o1n, vn3e1n;
-  vm3o1n.resize(3);
-  vm3e1n.resize(3);
-  vn3o1n.resize(3);
-  vn3e1n.resize(3);
-
-  std::vector<std::complex<double> > Ei, Hi, Es, Hs;
-  Ei.resize(3);
-  Hi.resize(3);
-  Es.resize(3);
-  Hs.resize(3);
-  for (i = 0; i < 3; i++) {
-    Ei[i] = std::complex<double>(0.0, 0.0);
-    Hi[i] = std::complex<double>(0.0, 0.0);
-    Es[i] = std::complex<double>(0.0, 0.0);
-    Hs[i] = std::complex<double>(0.0, 0.0);
-  }
-
-  std::vector<std::complex<double> > bj, by, bd;
-  bj.resize(nmax+1);
-  by.resize(nmax+1);
-  bd.resize(nmax+1);
-
-  // Calculate spherical Bessel and Hankel functions
-  sphericalBessel(Rho, nmax, bj, by, bd);
-
-  ci = std::complex<double>(0.0, 1.0);
-  for (n = 0; n < nmax; n++) {
-    n1 = n + 1;
-    rn = double(n + 1);
-
-    zn = bj[n1] + ci*by[n1];
-    xxip = Rho*(bj[n] + ci*by[n]) - rn*zn;
-
-    vm3o1n[0] = std::complex<double>(0.0, 0.0);
-    vm3o1n[1] = std::cos(Phi)*Pi[n]*zn;
-    vm3o1n[2] = -std::sin(Phi)*Tau[n]*zn;
-    vm3e1n[0] = std::complex<double>(0.0, 0.0);
-    vm3e1n[1] = -std::sin(Phi)*Pi[n]*zn;
-    vm3e1n[2] = -std::cos(Phi)*Tau[n]*zn;
-    vn3o1n[0] = std::sin(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
-    vn3o1n[1] = std::sin(Phi)*Tau[n]*xxip/Rho;
-    vn3o1n[2] = std::cos(Phi)*Pi[n]*xxip/Rho;
-    vn3e1n[0] = std::cos(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
-    vn3e1n[1] = std::cos(Phi)*Tau[n]*xxip/Rho;
-    vn3e1n[2] = -std::sin(Phi)*Pi[n]*xxip/Rho;
-
-    // scattered field: BH p.94 (4.45)
-    encap = std::pow(ci, rn)*(2.0*rn + 1.0)/(rn*rn + rn);
-    for (i = 0; i < 3; i++) {
-      Es[i] = Es[i] + encap*(ci*an[n]*vn3e1n[i] - bn[n]*vm3o1n[i]);
-      Hs[i] = Hs[i] + encap*(ci*bn[n]*vn3o1n[i] + an[n]*vm3e1n[i]);
-    }
-  }
-
-  // incident E field: BH p.89 (4.21); cf. p.92 (4.37), p.93 (4.38)
-  // basis unit vectors = er, etheta, ephi
-  std::complex<double> eifac = std::exp(std::complex<double>(0.0, Rho*std::cos(Theta)));
-
-  Ei[0] = eifac*std::sin(Theta)*std::cos(Phi);
-  Ei[1] = eifac*std::cos(Theta)*std::cos(Phi);
-  Ei[2] = -eifac*std::sin(Phi);
-
-  // magnetic field
-  double hffact = 1.0/(cc*mu);
-  for (i = 0; i < 3; i++) {
-    Hs[i] = hffact*Hs[i];
-  }
-
-  // incident H field: BH p.26 (2.43), p.89 (4.21)
-  std::complex<double> hffacta = hffact;
-  std::complex<double> hifac = eifac*hffacta;
-
-  Hi[0] = hifac*std::sin(Theta)*std::sin(Phi);
-  Hi[1] = hifac*std::cos(Theta)*std::sin(Phi);
-  Hi[2] = hifac*std::cos(Phi);
-
-  for (i = 0; i < 3; i++) {
-    // electric field E [V m-1] = EF*E0
-    E[i] = Ei[i] + Es[i];
-    H[i] = Hi[i] + Hs[i];
-  }
-}
-
-// Calculate an - equation (5)
-std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
-                             std::complex<double> PsiXL, std::complex<double> ZetaXL,
-                             std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
-
-  std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
-  std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
-
-  return Num/Denom;
-}
 
-// Calculate bn - equation (6)
-std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
-                             std::complex<double> PsiXL, std::complex<double> ZetaXL,
-                             std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
+  // ********************************************************************** //
+  // Modify refractive index of all layers                                  //
+  // ********************************************************************** //
+  void MultiLayerMie::SetLayersIndex(const std::vector< std::complex<double> >& index) {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    refr_index_ = index;
+  }
 
-  std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
-  std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
 
-  return Num/Denom;
-}
+  // ********************************************************************** //
+  // Modify coordinates for field calculation                               //
+  // ********************************************************************** //
+  void MultiLayerMie::SetFieldCoords(const std::vector< std::vector<double> >& coords) {
+    if (coords.size() != 3)
+      throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
+    if (coords[0].size() != coords[1].size() || coords[0].size() != coords[2].size())
+      throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
+    coords_ = coords;
+  }
 
-// Calculates S1 - equation (25a)
-std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
-                             double Pi, double Tau) {
 
-  return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
-}
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMie::SetPECLayer(int layer_position) {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    if (layer_position < 0)
+      throw std::invalid_argument("Error! Layers are numbered from 0!");
+    PEC_layer_position_ = layer_position;
+  }
 
-// Calculates S2 - equation (25b) (it's the same as (25a), just switches Pi and Tau)
-std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,
-                             double Pi, double Tau) {
 
-  return calc_S1(n, an, bn, Tau, Pi);
-}
+  // ********************************************************************** //
+  // Set maximun number of terms to be used                                 //
+  // ********************************************************************** //
+  void MultiLayerMie::SetMaxTerms(int nmax) {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    nmax_preset_ = nmax;
+  }
 
 
-//**********************************************************************************//
-// This function calculates the Riccati-Bessel functions (Psi and Zeta) for a       //
-// real argument (x).                                                               //
-// Equations (20a) - (21b)                                                          //
-//                                                                                  //
-// Input parameters:                                                                //
-//   x: Real argument to evaluate Psi and Zeta                                      //
-//   nmax: Maximum number of terms to calculate Psi and Zeta                        //
-//                                                                                  //
-// Output parameters:                                                               //
-//   Psi, Zeta: Riccati-Bessel functions                                            //
-//**********************************************************************************//
-void calcPsiZeta(double x, int nmax,
-                 std::vector<std::complex<double> > D1,
-                 std::vector<std::complex<double> > D3,
-                 std::vector<std::complex<double> >& Psi,
-                 std::vector<std::complex<double> >& Zeta) {
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double MultiLayerMie::GetSizeParameter() {
+    if (size_param_.size() > 0)
+      return size_param_.back();
+    else
+      return 0;
+  }
 
-  int n;
 
-  //Upward recurrence for Psi and Zeta - equations (20a) - (21b)
-  Psi[0] = std::complex<double>(sin(x), 0);
-  Zeta[0] = std::complex<double>(sin(x), -cos(x));
-  for (n = 1; n <= nmax; n++) {
-    Psi[n] = Psi[n - 1]*(n/x - D1[n - 1]);
-    Zeta[n] = Zeta[n - 1]*(n/x - D3[n - 1]);
+  // ********************************************************************** //
+  // Clear layer information                                                //
+  // ********************************************************************** //
+  void MultiLayerMie::ClearLayers() {
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+    size_param_.clear();
+    refr_index_.clear();
   }
-}
-
-//**********************************************************************************//
-// This function calculates the logarithmic derivatives of the Riccati-Bessel       //
-// functions (D1 and D3) for a complex argument (z).                                //
-// Equations (16a), (16b) and (18a) - (18d)                                         //
-//                                                                                  //
-// Input parameters:                                                                //
-//   z: Complex argument to evaluate D1 and D3                                      //
-//   nmax: Maximum number of terms to calculate D1 and D3                           //
-//                                                                                  //
-// Output parameters:                                                               //
-//   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions                //
-//**********************************************************************************//
-void calcD1D3(std::complex<double> z, int nmax,
-              std::vector<std::complex<double> >& D1,
-              std::vector<std::complex<double> >& D3) {
 
-  int n;
-  std::complex<double> nz, PsiZeta;
 
-  // Downward recurrence for D1 - equations (16a) and (16b)
-  D1[nmax] = std::complex<double>(0.0, 0.0);
-  for (n = nmax; n > 0; n--) {
-    nz = double(n)/z;
-    D1[n - 1] = nz - 1.0/(D1[n] + nz);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  //                         Computational core
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+
+
+  // ********************************************************************** //
+  // Calculate calcNstop - equation (17)                                        //
+  // ********************************************************************** //
+  void MultiLayerMie::calcNstop() {
+    const double& xL = size_param_.back();
+    if (xL <= 8) {
+      nmax_ = round(xL + 4.0*pow(xL, 1.0/3.0) + 1);
+    } else if (xL <= 4200) {
+      nmax_ = round(xL + 4.05*pow(xL, 1.0/3.0) + 2);
+    } else {
+      nmax_ = round(xL + 4.0*pow(xL, 1.0/3.0) + 2);
+    }
   }
 
-  // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
-  PsiZeta = 0.5*(1.0 - std::complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))*exp(-2.0*z.imag()));
-  D3[0] = std::complex<double>(0.0, 1.0);
-  for (n = 1; n <= nmax; n++) {
-    nz = double(n)/z;
-    PsiZeta = PsiZeta*(nz - D1[n - 1])*(nz - D3[n - 1]);
-    D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta;
-  }
-}
 
-//**********************************************************************************//
-// This function calculates Pi and Tau for all values of Theta.                     //
-// Equations (26a) - (26c)                                                          //
-//                                                                                  //
-// Input parameters:                                                                //
-//   nmax: Maximum number of terms to calculate Pi and Tau                          //
-//   nTheta: Number of scattering angles                                            //
-//   Theta: Array containing all the scattering angles where the scattering         //
-//          amplitudes will be calculated                                           //
-//                                                                                  //
-// Output parameters:                                                               //
-//   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
-//**********************************************************************************//
-void calcPiTau(int nmax, double Theta, std::vector<double>& Pi, std::vector<double>& Tau) {
-
-  int n;
-  //****************************************************//
-  // Equations (26a) - (26c)                            //
-  //****************************************************//
-  // Initialize Pi and Tau
-  Pi[0] = 1.0;
-  Tau[0] = cos(Theta);
-  // Calculate the actual values
-  if (nmax > 1) {
-    Pi[1] = 3*Tau[0]*Pi[0];
-    Tau[1] = 2*Tau[0]*Pi[1] - 3*Pi[0];
-    for (n = 2; n < nmax; n++) {
-      Pi[n] = ((n + n + 1)*Tau[0]*Pi[n - 1] - (n + 1)*Pi[n - 2])/n;
-      Tau[n] = (n + 1)*Tau[0]*Pi[n] - (n + 2)*Pi[n - 1];
-    }
-  }
-}
+  // ********************************************************************** //
+  // Maximum number of terms required for the calculation                   //
+  // ********************************************************************** //
+  void MultiLayerMie::calcNmax(int first_layer) {
+    int ri, riM1;
+    const std::vector<double>& x = size_param_;
+    const std::vector<std::complex<double> >& m = refr_index_;
+    calcNstop();  // Set initial nmax_ value
+    for (int i = first_layer; i < x.size(); i++) {
+      if (i > PEC_layer_position_)
+        ri = round(std::abs(x[i]*m[i]));
+      else
+        ri = 0;
+      nmax_ = std::max(nmax_, ri);
+      // first layer is pec, if pec is present
+      if ((i > first_layer) && ((i - 1) > PEC_layer_position_))
+        riM1 = round(std::abs(x[i - 1]* m[i]));
+      else
+        riM1 = 0;
+      nmax_ = std::max(nmax_, riM1);
+    }
+    nmax_ += 15;  // Final nmax_ value
+  }
 
-//**********************************************************************************//
-// This function calculates 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(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
-                std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
-  //************************************************************************//
-  // Calculate the index of the first layer. It can be either 0 (default)   //
-  // or the index of the outermost PEC layer. In the latter case all layers //
-  // below the PEC are discarded.                                           //
-  //************************************************************************//
 
-  int fl = (pl > 0) ? pl : 0;
+  //**********************************************************************************//
+  // This function calculates the spherical Bessel (jn) and Hankel (h1n) functions    //
+  // and their derivatives for a given complex value z. See pag. 87 B&H.              //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   z: Complex argument to evaluate jn and h1n                                     //
+  //   nmax_: Maximum number of terms to calculate jn and h1n                         //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   jn, h1n: Spherical Bessel and Hankel functions                                 //
+  //   jnp, h1np: Derivatives of the spherical Bessel and Hankel functions            //
+  //                                                                                  //
+  // The implementation follows the algorithm by I.J. Thompson and A.R. Barnett,      //
+  // Comp. Phys. Comm. 47 (1987) 245-257.                                             //
+  //                                                                                  //
+  // Complex spherical Bessel functions from n=0..nmax_-1 for z in the upper half     //
+  // plane (Im(z) > -3).                                                              //
+  //                                                                                  //
+  //     j[n]   = j/n(z)                Regular solution: j[0]=sin(z)/z               //
+  //     j'[n]  = d[j/n(z)]/dz                                                        //
+  //     h1[n]  = h[0]/n(z)             Irregular Hankel function:                    //
+  //     h1'[n] = d[h[0]/n(z)]/dz                h1[0] = j0(z) + i*y0(z)              //
+  //                                                   = (sin(z)-i*cos(z))/z          //
+  //                                                   = -i*exp(i*z)/z                //
+  // Using complex CF1, and trigonometric forms for n=0 solutions.                    //
+  //**********************************************************************************//
+  void MultiLayerMie::sbesjh(std::complex<double> z,
+                             std::vector<std::complex<double> >& jn, std::vector<std::complex<double> >& jnp,
+                             std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np) {
+    const int limit = 20000;
+    const double accur = 1.0e-12;
+    const double tm30 = 1e-30;
+    const std::complex<double> c_i(0.0, 1.0);
+
+    double absc;
+    std::complex<double> zi, w;
+    std::complex<double> pl, f, b, d, c, del, jn0;
+
+    absc = std::abs(std::real(z)) + std::abs(std::imag(z));
+    if ((absc < accur) || (std::imag(z) < -3.0)) {
+      throw std::invalid_argument("TODO add error description for condition if ((absc < accur) || (std::imag(z) < -3.0))");
+    }
 
-  if (nmax <= 0) {
-    nmax = Nmax(L, fl, pl, x, m);
-  }
+    zi = 1.0/z;
+    w = zi + zi;
 
-  std::complex<double> z1, z2;
-  std::complex<double> Num, Denom;
-  std::complex<double> G1, G2;
-  std::complex<double> Temp;
+    pl = double(nmax_)*zi;
 
-  int n, l;
+    f = pl + zi;
+    b = f + f + zi;
+    d = 0.0;
+    c = f;
+    for (int l = 0; l < limit; l++) {
+      d = b - d;
+      c = b - 1.0/c;
 
-  //**************************************************************************//
-  // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which  //
-  // means that index = layer number - 1 or index = n - 1. The only exception //
-  // are the arrays for representing D1, D3 and Q because they need a value   //
-  // for the index 0 (zero), hence it is important to consider this shift     //
-  // between different arrays. The change was done to optimize memory usage.  //
-  //**************************************************************************//
+      absc = std::abs(std::real(d)) + std::abs(std::imag(d));
+      if (absc < tm30) {
+        d = tm30;
+      }
 
-  // Allocate memory to the arrays
-  std::vector<std::vector<std::complex<double> > > D1_mlxl, D1_mlxlM1;
-  D1_mlxl.resize(L);
-  D1_mlxlM1.resize(L);
+      absc = std::abs(std::real(c)) + std::abs(std::imag(c));
+      if (absc < tm30) {
+        c = tm30;
+      }
 
-  std::vector<std::vector<std::complex<double> > > D3_mlxl, D3_mlxlM1;
-  D3_mlxl.resize(L);
-  D3_mlxlM1.resize(L);
+      d = 1.0/d;
+      del = d*c;
+      f = f*del;
+      b += w;
 
-  std::vector<std::vector<std::complex<double> > > Q;
-  Q.resize(L);
+      absc = std::abs(std::real(del - 1.0)) + std::abs(std::imag(del - 1.0));
 
-  std::vector<std::vector<std::complex<double> > > Ha, Hb;
-  Ha.resize(L);
-  Hb.resize(L);
+      // We have obtained the desired accuracy
+      if (absc < accur) break;
+    }
 
-  for (l = 0; l < L; l++) {
-    D1_mlxl[l].resize(nmax + 1);
-    D1_mlxlM1[l].resize(nmax + 1);
+    if (absc > accur) {
+      throw std::invalid_argument("We were not able to obtain the desired accuracy (MultiLayerMie::sbesjh)");
+    }
 
-    D3_mlxl[l].resize(nmax + 1);
-    D3_mlxlM1[l].resize(nmax + 1);
+    jn[nmax_ - 1] = tm30;
+    jnp[nmax_ - 1] = f*jn[nmax_ - 1];
 
-    Q[l].resize(nmax + 1);
+    // Downward recursion to n=0 (N.B.  Coulomb Functions)
+    for (int n = nmax_ - 2; n >= 0; n--) {
+      jn[n] = pl*jn[n + 1] + jnp[n + 1];
+      jnp[n] = pl*jn[n] - jn[n + 1];
+      pl = pl - zi;
+    }
 
-    Ha[l].resize(nmax);
-    Hb[l].resize(nmax);
+    // Calculate the n=0 Bessel Functions
+    jn0 = zi*std::sin(z);
+    h1n[0] = -c_i*zi*std::exp(c_i*z);
+    h1np[0] = h1n[0]*(c_i - zi);
+
+    // Rescale j[n], j'[n], converting to spherical Bessel functions.
+    // Recur   h1[n], h1'[n] as spherical Bessel functions.
+    w = 1.0/jn[0];
+    pl = zi;
+    for (int n = 0; n < nmax_; n++) {
+      jn[n] = jn0*w*jn[n];
+      jnp[n] = jn0*w*jnp[n] - zi*jn[n];
+      if (n > 0) {
+        h1n[n] = (pl - zi)*h1n[n - 1] - h1np[n - 1];
+
+        // check if hankel is increasing (upward stable)
+        if (std::abs(h1n[n]) < std::abs(h1n[n - 1])) {
+          throw std::invalid_argument("Error: Hankel not increasing! (MultiLayerMie::sbesjh)");
+        }
+
+        pl += zi;
+
+        h1np[n] = -pl*h1n[n] + h1n[n - 1];
+      }
+    }
   }
 
-  an.resize(nmax);
-  bn.resize(nmax);
+  // ********************************************************************** //
+  // Calculate an - equation (5)                                            //
+  // ********************************************************************** //
+  std::complex<double> MultiLayerMie::calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
+                                              std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                                              std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
 
-  std::vector<std::complex<double> > D1XL, D3XL;
-  D1XL.resize(nmax + 1);
-  D3XL.resize(nmax + 1);
+    std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
+    std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
 
+    return Num/Denom;
+  }
 
-  std::vector<std::complex<double> > PsiXL, ZetaXL;
-  PsiXL.resize(nmax + 1);
-  ZetaXL.resize(nmax + 1);
 
-  //*************************************************//
-  // Calculate D1 and D3 for z1 in the first layer   //
-  //*************************************************//
-  if (fl == pl) {  // PEC layer
-    for (n = 0; n <= nmax; n++) {
-      D1_mlxl[fl][n] = std::complex<double>(0.0, -1.0);
-      D3_mlxl[fl][n] = std::complex<double>(0.0, 1.0);
-    }
-  } else { // Regular layer
-    z1 = x[fl]* m[fl];
+  // ********************************************************************** //
+  // Calculate bn - equation (6)                                            //
+  // ********************************************************************** //
+  std::complex<double> MultiLayerMie::calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
+                                              std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                                              std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
+
+    std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
+    std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
 
-    // Calculate D1 and D3
-    calcD1D3(z1, nmax, D1_mlxl[fl], D3_mlxl[fl]);
+    return Num/Denom;
   }
 
-  //******************************************************************//
-  // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
-  //******************************************************************//
-  for (n = 0; n < nmax; n++) {
-    Ha[fl][n] = D1_mlxl[fl][n + 1];
-    Hb[fl][n] = D1_mlxl[fl][n + 1];
+
+  // ********************************************************************** //
+  // Calculates S1 - equation (25a)                                         //
+  // ********************************************************************** //
+  std::complex<double> MultiLayerMie::calc_S1(int n, std::complex<double> an, std::complex<double> bn,
+                                              double Pi, double Tau) {
+    return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
   }
 
-  //*****************************************************//
-  // Iteration from the second layer to the last one (L) //
-  //*****************************************************//
-  for (l = fl + 1; l < L; l++) {
-    //************************************************************//
-    //Calculate D1 and D3 for z1 and z2 in the layers fl+1..L     //
-    //************************************************************//
-    z1 = x[l]*m[l];
-    z2 = x[l - 1]*m[l];
 
-    //Calculate D1 and D3 for z1
-    calcD1D3(z1, nmax, D1_mlxl[l], D3_mlxl[l]);
+  // ********************************************************************** //
+  // Calculates S2 - equation (25b) (it's the same as (25a), just switches  //
+  // Pi and Tau)                                                            //
+  // ********************************************************************** //
+  std::complex<double> MultiLayerMie::calc_S2(int n, std::complex<double> an, std::complex<double> bn,
+                                              double Pi, double Tau) {
+    return calc_S1(n, an, bn, Tau, Pi);
+  }
 
-    //Calculate D1 and D3 for z2
-    calcD1D3(z2, nmax, D1_mlxlM1[l], D3_mlxlM1[l]);
 
-    //*********************************************//
-    //Calculate Q, Ha and Hb in the layers fl+1..L //
-    //*********************************************//
+  //**********************************************************************************//
+  // This function calculates the Riccati-Bessel functions (Psi and Zeta) for a       //
+  // real argument (x).                                                               //
+  // Equations (20a) - (21b)                                                          //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   x: Real argument to evaluate Psi and Zeta                                      //
+  //   nmax: Maximum number of terms to calculate Psi and Zeta                        //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   Psi, Zeta: Riccati-Bessel functions                                            //
+  //**********************************************************************************//
+  void MultiLayerMie::calcPsiZeta(std::complex<double> z,
+                                  std::vector<std::complex<double> > D1,
+                                  std::vector<std::complex<double> > D3,
+                                  std::vector<std::complex<double> >& Psi,
+                                  std::vector<std::complex<double> >& Zeta) {
+
+    //Upward recurrence for Psi and Zeta - equations (20a) - (21b)
+    std::complex<double> c_i(0.0, 1.0);
+    Psi[0] = std::sin(z);
+    Zeta[0] = std::sin(z) - c_i*std::cos(z);
+    for (int n = 1; n <= nmax_; n++) {
+      Psi[n] = Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
+      Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
+    }
 
-    // Upward recurrence for Q - equations (19a) and (19b)
-    Num = exp(-2.0*(z1.imag() - z2.imag()))*std::complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
-    Denom = std::complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
-    Q[l][0] = Num/Denom;
+  }
 
-    for (n = 1; n <= nmax; n++) {
-      Num = (z1*D1_mlxl[l][n] + double(n))*(double(n) - z1*D3_mlxl[l][n - 1]);
-      Denom = (z2*D1_mlxlM1[l][n] + double(n))*(double(n) - z2*D3_mlxlM1[l][n - 1]);
 
-      Q[l][n] = (((x[l - 1]*x[l - 1])/(x[l]*x[l])* Q[l][n - 1])*Num)/Denom;
+  //**********************************************************************************//
+  // This function calculates the logarithmic derivatives of the Riccati-Bessel       //
+  // functions (D1 and D3) for a complex argument (z).                                //
+  // Equations (16a), (16b) and (18a) - (18d)                                         //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   z: Complex argument to evaluate D1 and D3                                      //
+  //   nmax_: Maximum number of terms to calculate D1 and D3                          //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions                //
+  //**********************************************************************************//
+  void MultiLayerMie::calcD1D3(const std::complex<double> z,
+                               std::vector<std::complex<double> >& D1,
+                               std::vector<std::complex<double> >& D3) {
+
+    // Downward recurrence for D1 - equations (16a) and (16b)
+    D1[nmax_] = std::complex<double>(0.0, 0.0);
+    const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
+
+    for (int n = nmax_; n > 0; n--) {
+      D1[n - 1] = double(n)*zinv - 1.0/(D1[n] + double(n)*zinv);
     }
 
-    // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
-    for (n = 1; n <= nmax; n++) {
-      //Ha
-      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
-        G1 = -D1_mlxlM1[l][n];
-        G2 = -D3_mlxlM1[l][n];
-      } else {
-        G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[l][n]);
-        G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[l][n]);
+    if (std::abs(D1[0]) > 100000.0)
+      throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
+
+    // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
+    PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
+                       *std::exp(-2.0*z.imag()));
+    D3[0] = std::complex<double>(0.0, 1.0);
+    for (int n = 1; n <= nmax_; n++) {
+      PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
+        *(static_cast<double>(n)*zinv- D3[n - 1]);
+      D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
+    }
+  }
+
+
+  //**********************************************************************************//
+  // This function calculates Pi and Tau for a given value of cos(Theta).             //
+  // Equations (26a) - (26c)                                                          //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   nmax_: Maximum number of terms to calculate Pi and Tau                         //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
+  //**********************************************************************************//
+  void MultiLayerMie::calcPiTau(const double& costheta,
+                                std::vector<double>& Pi, std::vector<double>& Tau) {
+
+    int n;
+    //****************************************************//
+    // Equations (26a) - (26c)                            //
+    //****************************************************//
+    // Initialize Pi and Tau
+    Pi[0] = 1.0;
+    Tau[0] = costheta;
+    // Calculate the actual values
+    if (nmax_ > 1) {
+      Pi[1] = 3*costheta*Pi[0];
+      Tau[1] = 2*costheta*Pi[1] - 3*Pi[0];
+      for (n = 2; n < nmax_; n++) {
+        Pi[n] = ((n + n + 1)*costheta*Pi[n - 1] - (n + 1)*Pi[n - 2])/n;
+        Tau[n] = (n + 1)*costheta*Pi[n] - (n + 2)*Pi[n - 1];
       }
+    }
+  }  // end of MultiLayerMie::calcPiTau(...)
+
+  void MultiLayerMie::calcSpherHarm(const double Rho, const double Phi, const double Theta,
+                                    const std::complex<double>& zn, const std::complex<double>& dzn,
+                                    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) {
+
+    // using eq 4.50 in BH
+    std::complex<double> c_zero(0.0, 0.0);
+    std::complex<double> deriv = Rho*dzn + zn;
+
+    using std::sin;
+    using std::cos;
+    Mo1n[0] = c_zero;
+    Mo1n[1] = cos(Phi)*Pi*zn;
+    Mo1n[2] = -sin(Phi)*Tau*zn;
+    Me1n[0] = c_zero;
+    Me1n[1] = -sin(Phi)*Pi*zn;
+    Me1n[2] = -cos(Phi)*Tau*zn;
+    No1n[0] = sin(Phi)*(n*n + n)*sin(Theta)*Pi*zn/Rho;
+    No1n[1] = sin(Phi)*Tau*deriv/Rho;
+    No1n[2] = cos(Phi)*Pi*deriv/Rho;
+    Ne1n[0] = cos(Phi)*(n*n + n)*sin(Theta)*Pi*zn/Rho;
+    Ne1n[1] = cos(Phi)*Tau*deriv/Rho;
+    Ne1n[2] = -sin(Phi)*Pi*deriv/Rho;
+  }  // end of MultiLayerMie::calcSpherHarm(...)
+
+
+  //**********************************************************************************//
+  // This function calculates 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                 //
+  //**********************************************************************************//
+  void MultiLayerMie::ExtScattCoeffs() {
+
+    areExtCoeffsCalc_ = false;
+
+    const std::vector<double>& x = size_param_;
+    const std::vector<std::complex<double> >& m = refr_index_;
+    const int& pl = PEC_layer_position_;
+    const int L = refr_index_.size();
+
+    //************************************************************************//
+    // Calculate the index of the first layer. It can be either 0 (default)   //
+    // or the index of the outermost PEC layer. In the latter case all layers //
+    // below the PEC are discarded.                                           //
+    // ***********************************************************************//
+    // TODO, is it possible for PEC to have a zero index? If yes than
+    // is should be:
+    // int fl = (pl > - 1) ? pl : 0;
+    // This will give the same result, however, it corresponds the
+    // logic - if there is PEC, than first layer is PEC.
+    // Well, I followed the logic: First layer is always zero unless it has
+    // an upper PEC layer.
+    int fl = (pl > 0) ? pl : 0;
+    if (nmax_preset_ <= 0) calcNmax(fl);
+    else nmax_ = nmax_preset_;
+
+    std::complex<double> z1, z2;
+    //**************************************************************************//
+    // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which  //
+    // means that index = layer number - 1 or index = n - 1. The only exception //
+    // are the arrays for representing D1, D3 and Q because they need a value   //
+    // for the index 0 (zero), hence it is important to consider this shift     //
+    // between different arrays. The change was done to optimize memory usage.  //
+    //**************************************************************************//
+    // Allocate memory to the arrays
+    std::vector<std::complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
+                                       D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
+
+    std::vector<std::vector<std::complex<double> > > Q(L), Ha(L), Hb(L);
+
+    for (int l = 0; l < L; l++) {
+      Q[l].resize(nmax_ + 1);
+      Ha[l].resize(nmax_);
+      Hb[l].resize(nmax_);
+    }
 
-      Temp = Q[l][n]*G1;
+    an_.resize(nmax_);
+    bn_.resize(nmax_);
+    PsiZeta_.resize(nmax_ + 1);
 
-      Num = (G2*D1_mlxl[l][n]) - (Temp*D3_mlxl[l][n]);
-      Denom = G2 - Temp;
+    std::vector<std::complex<double> > D1XL(nmax_ + 1), D3XL(nmax_ + 1),
+                                       PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
 
-      Ha[l][n - 1] = Num/Denom;
+    //*************************************************//
+    // Calculate D1 and D3 for z1 in the first layer   //
+    //*************************************************//
+    if (fl == pl) {  // PEC layer
+      for (int n = 0; n <= nmax_; n++) {
+        D1_mlxl[n] = std::complex<double>(0.0, - 1.0);
+        D3_mlxl[n] = std::complex<double>(0.0, 1.0);
+      }
+    } else { // Regular layer
+      z1 = x[fl]* m[fl];
+      // Calculate D1 and D3
+      calcD1D3(z1, D1_mlxl, D3_mlxl);
+    }
 
-      //Hb
-      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
-        G1 = Hb[l - 1][n - 1];
-        G2 = Hb[l - 1][n - 1];
+    //******************************************************************//
+    // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
+    //******************************************************************//
+    for (int n = 0; n < nmax_; n++) {
+      Ha[fl][n] = D1_mlxl[n + 1];
+      Hb[fl][n] = D1_mlxl[n + 1];
+    }
+    //*****************************************************//
+    // Iteration from the second layer to the last one (L) //
+    //*****************************************************//
+    std::complex<double> Temp, Num, Denom;
+    std::complex<double> G1, G2;
+    for (int l = fl + 1; l < L; l++) {
+      //************************************************************//
+      //Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L     //
+      //************************************************************//
+      z1 = x[l]*m[l];
+      z2 = x[l - 1]*m[l];
+      //Calculate D1 and D3 for z1
+      calcD1D3(z1, D1_mlxl, D3_mlxl);
+      //Calculate D1 and D3 for z2
+      calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
+
+      //*********************************************//
+      //Calculate Q, Ha and Hb in the layers fl + 1..L //
+      //*********************************************//
+      // Upward recurrence for Q - equations (19a) and (19b)
+      Num = std::exp(-2.0*(z1.imag() - z2.imag()))
+           *std::complex<double>(std::cos(-2.0*z2.real()) - std::exp(-2.0*z2.imag()), std::sin(-2.0*z2.real()));
+      Denom = std::complex<double>(std::cos(-2.0*z1.real()) - std::exp(-2.0*z1.imag()), std::sin(-2.0*z1.real()));
+      Q[l][0] = Num/Denom;
+      for (int n = 1; n <= nmax_; n++) {
+        Num = (z1*D1_mlxl[n] + double(n))*(double(n) - z1*D3_mlxl[n - 1]);
+        Denom = (z2*D1_mlxlM1[n] + double(n))*(double(n) - z2*D3_mlxlM1[n - 1]);
+        Q[l][n] = ((pow2(x[l - 1]/x[l])* Q[l][n - 1])*Num)/Denom;
+      }
+      // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
+      for (int n = 1; n <= nmax_; n++) {
+        //Ha
+        if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+          G1 = -D1_mlxlM1[n];
+          G2 = -D3_mlxlM1[n];
+        } else {
+          G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[n]);
+          G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[n]);
+        }  // end of if PEC
+        Temp = Q[l][n]*G1;
+        Num = (G2*D1_mlxl[n]) - (Temp*D3_mlxl[n]);
+        Denom = G2 - Temp;
+        Ha[l][n - 1] = Num/Denom;
+        //Hb
+        if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+          G1 = Hb[l - 1][n - 1];
+          G2 = Hb[l - 1][n - 1];
+        } else {
+          G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[n]);
+          G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[n]);
+        }  // end of if PEC
+
+        Temp = Q[l][n]*G1;
+        Num = (G2*D1_mlxl[n]) - (Temp* D3_mlxl[n]);
+        Denom = (G2- Temp);
+        Hb[l][n - 1] = (Num/ Denom);
+      }  // end of for Ha and Hb terms
+    }  // end of for layers iteration
+
+    //**************************************//
+    //Calculate D1, D3, Psi and Zeta for XL //
+    //**************************************//
+    // Calculate D1XL and D3XL
+    calcD1D3(x[L - 1], D1XL, D3XL);
+
+    // Calculate PsiXL and ZetaXL
+    calcPsiZeta(x[L - 1], D1XL, D3XL, PsiXL, ZetaXL);
+    //*********************************************************************//
+    // Finally, we calculate the scattering coefficients (an and bn) and   //
+    // the angular functions (Pi and Tau). Note that for these arrays the  //
+    // first layer is 0 (zero), in future versions all arrays will follow  //
+    // this convention to save memory. (13 Nov, 2014)                      //
+    //*********************************************************************//
+    for (int n = 0; n < nmax_; n++) {
+      //********************************************************************//
+      //Expressions for calculating an and bn coefficients are not valid if //
+      //there is only one PEC layer (ie, for a simple PEC sphere).          //
+      //********************************************************************//
+      if (pl < (L - 1)) {
+        an_[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+        bn_[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
       } else {
-        G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[l][n]);
-        G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[l][n]);
+        an_[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+        bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
+      }
+    }  // end of for an and bn terms
+    areExtCoeffsCalc_ = true;
+  }  // end of void MultiLayerMie::ExtScattCoeffs(...)
+
+
+  //**********************************************************************************//
+  // This function calculates the actual scattering parameters and amplitudes         //
+  //                                                                                  //
+  // 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]     //
+  //   nTheta: Number of scattering angles                                            //
+  //   Theta: Array containing all the scattering angles where the scattering         //
+  //          amplitudes will be calculated                                           //
+  //   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:                                                               //
+  //   Qext: Efficiency factor for extinction                                         //
+  //   Qsca: Efficiency factor for scattering                                         //
+  //   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
+  //   Qbk: Efficiency factor for backscattering                                      //
+  //   Qpr: Efficiency factor for the radiation pressure                              //
+  //   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
+  //   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
+  //   S1, S2: Complex scattering amplitudes                                          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  void MultiLayerMie::RunMieCalculation() {
+    if (size_param_.size() != refr_index_.size())
+      throw std::invalid_argument("Each size parameter should have only one index!");
+    if (size_param_.size() == 0)
+      throw std::invalid_argument("Initialize model first!");
+
+    const std::vector<double>& x = size_param_;
+
+    areIntCoeffsCalc_ = false;
+    areExtCoeffsCalc_ = false;
+    isMieCalculated_ = false;
+
+    // Calculate scattering coefficients
+    ExtScattCoeffs();
+
+//    for (int i = 0; i < nmax_; i++) {
+//      printf("a[%i] = %g, %g; b[%i] = %g, %g\n", i, an_[i].real(), an_[i].imag(), i, bn_[i].real(), bn_[i].imag());
+//    }
+
+    if (!areExtCoeffsCalc_)
+      throw std::invalid_argument("Calculation of scattering coefficients failed!");
+
+    // Initialize the scattering parameters
+    Qext_ = 0;
+    Qsca_ = 0;
+    Qabs_ = 0;
+    Qbk_ = 0;
+    Qpr_ = 0;
+    asymmetry_factor_ = 0;
+    albedo_ = 0;
+    Qsca_ch_.clear();
+    Qext_ch_.clear();
+    Qabs_ch_.clear();
+    Qbk_ch_.clear();
+    Qpr_ch_.clear();
+    Qsca_ch_.resize(nmax_ - 1);
+    Qext_ch_.resize(nmax_ - 1);
+    Qabs_ch_.resize(nmax_ - 1);
+    Qbk_ch_.resize(nmax_ - 1);
+    Qpr_ch_.resize(nmax_ - 1);
+    Qsca_ch_norm_.resize(nmax_ - 1);
+    Qext_ch_norm_.resize(nmax_ - 1);
+    Qabs_ch_norm_.resize(nmax_ - 1);
+    Qbk_ch_norm_.resize(nmax_ - 1);
+    Qpr_ch_norm_.resize(nmax_ - 1);
+
+    // Initialize the scattering amplitudes
+    std::vector<std::complex<double> > tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
+    S1_.swap(tmp1);
+    S2_ = S1_;
+
+    std::vector<double> Pi(nmax_), Tau(nmax_);
+
+    std::complex<double> Qbktmp(0.0, 0.0);
+    std::vector< std::complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
+    // By using downward recurrence we avoid loss of precision due to float rounding errors
+    // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
+    //      http://en.wikipedia.org/wiki/Loss_of_significance
+    for (int i = nmax_ - 2; i >= 0; i--) {
+      const int n = i + 1;
+      // Equation (27)
+      Qext_ch_norm_[i] = (an_[i].real() + bn_[i].real());
+      Qext_ch_[i] = (n + n + 1.0)*Qext_ch_norm_[i];
+      //Qext_ch_[i] = (n + n + 1)*(an_[i].real() + bn_[i].real());
+      Qext_ += Qext_ch_[i];
+      // Equation (28)
+      Qsca_ch_norm_[i] = (an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
+                          + bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
+      Qsca_ch_[i] = (n + n + 1.0)*Qsca_ch_norm_[i];
+      Qsca_ += Qsca_ch_[i];
+      // Qsca_ch_[i] += (n + n + 1)*(an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
+      //                             + bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
+
+      // Equation (29) TODO We must check carefully this equation. If we
+      // remove the typecast to double then the result changes. Which is
+      // the correct one??? Ovidio (2014/12/10) With cast ratio will
+      // give double, without cast (n + n + 1)/(n*(n + 1)) will be
+      // rounded to integer. Tig (2015/02/24)
+      Qpr_ch_[i]=((n*(n + 2)/(n + 1))*((an_[i]*std::conj(an_[n]) + bn_[i]*std::conj(bn_[n])).real())
+               + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*std::conj(bn_[i])).real());
+      Qpr_ += Qpr_ch_[i];
+      // Equation (33)
+      Qbktmp_ch[i] = (double)(n + n + 1)*(1 - 2*(n % 2))*(an_[i]- bn_[i]);
+      Qbktmp += Qbktmp_ch[i];
+      // Calculate the scattering amplitudes (S1 and S2)    //
+      // Equations (25a) - (25b)                            //
+      for (int t = 0; t < theta_.size(); t++) {
+        calcPiTau(std::cos(theta_[t]), Pi, Tau);
+
+        S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
+        S2_[t] += calc_S2(n, an_[i], bn_[i], Pi[i], Tau[i]);
       }
+    }
+    double x2 = pow2(x.back());
+    Qext_ = 2.0*(Qext_)/x2;                                 // Equation (27)
+    for (double& Q : Qext_ch_) Q = 2.0*Q/x2;
+    Qsca_ = 2.0*(Qsca_)/x2;                                 // Equation (28)
+    for (double& Q : Qsca_ch_) Q = 2.0*Q/x2;
+    //for (double& Q : Qsca_ch_norm_) Q = 2.0*Q/x2;
+    Qpr_ = Qext_ - 4.0*(Qpr_)/x2;                           // Equation (29)
+    for (int i = 0; i < nmax_ - 1; ++i) Qpr_ch_[i] = Qext_ch_[i] - 4.0*Qpr_ch_[i]/x2;
+
+    Qabs_ = Qext_ - Qsca_;                                // Equation (30)
+    for (int i = 0; i < nmax_ - 1; ++i) {
+      Qabs_ch_[i] = Qext_ch_[i] - Qsca_ch_[i];
+      Qabs_ch_norm_[i] = Qext_ch_norm_[i] - Qsca_ch_norm_[i];
+    }
 
-      Temp = Q[l][n]*G1;
+    albedo_ = Qsca_/Qext_;                              // Equation (31)
+    asymmetry_factor_ = (Qext_ - Qpr_)/Qsca_;                          // Equation (32)
 
-      Num = (G2*D1_mlxl[l][n]) - (Temp* D3_mlxl[l][n]);
-      Denom = (G2- Temp);
+    Qbk_ = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2;    // Equation (33)
 
-      Hb[l][n - 1] = (Num/ Denom);
-    }
+    isMieCalculated_ = true;
   }
 
-  //**************************************//
-  //Calculate D1, D3, Psi and Zeta for XL //
-  //**************************************//
 
-  // Calculate D1XL and D3XL
-  calcD1D3(x[L - 1], nmax, D1XL, D3XL);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMie::IntScattCoeffs() {
+    if (!areExtCoeffsCalc_)
+      throw std::invalid_argument("(IntScattCoeffs) You should calculate external coefficients first!");
+
+    areIntCoeffsCalc_ = false;
+
+    std::complex<double> c_one(1.0, 0.0);
+    std::complex<double> c_zero(0.0, 0.0);
+
+    const int L = refr_index_.size();
+
+    // we need to fill
+    // std::vector< std::vector<std::complex<double> > > anl_, bnl_, cnl_, dnl_;
+    //     for n = [0..nmax_) and for l=[L..0)
+    // TODO: to decrease cache miss outer loop is with n and inner with reversed l
+    // at the moment outer is forward l and inner in n
+    anl_.resize(L + 1);
+    bnl_.resize(L + 1);
+    cnl_.resize(L + 1);
+    dnl_.resize(L + 1);
+    for (auto& element:anl_) element.resize(nmax_);
+    for (auto& element:bnl_) element.resize(nmax_);
+    for (auto& element:cnl_) element.resize(nmax_);
+    for (auto& element:dnl_) element.resize(nmax_);
+
+    // Yang, paragraph under eq. A3
+    // a^(L + 1)_n = a_n, d^(L + 1) = 1 ...
+    for (int i = 0; i < nmax_; ++i) {
+      anl_[L][i] = an_[i];
+      bnl_[L][i] = bn_[i];
+      cnl_[L][i] = c_one;
+      dnl_[L][i] = c_one;
+    }
 
-  // Calculate PsiXL and ZetaXL
-  calcPsiZeta(x[L - 1], nmax, D1XL, D3XL, PsiXL, ZetaXL);
+    std::vector<std::complex<double> > z(L), z1(L);
+    for (int i = 0; i < L - 1; ++i) {
+      z[i] = size_param_[i]*refr_index_[i];
+      z1[i] = size_param_[i]*refr_index_[i + 1];
+    }
+    z[L - 1] = size_param_[L - 1]*refr_index_[L - 1];
+    z1[L - 1] = size_param_[L - 1];
+    std::vector< std::vector<std::complex<double> > > D1z(L), D1z1(L), D3z(L), D3z1(L);
+    std::vector< std::vector<std::complex<double> > > Psiz(L), Psiz1(L), Zetaz(L), Zetaz1(L);
+    for (int l = 0; l < L; ++l) {
+      D1z[l].resize(nmax_ + 1);
+      D1z1[l].resize(nmax_ + 1);
+      D3z[l].resize(nmax_ + 1);
+      D3z1[l].resize(nmax_ + 1);
+      Psiz[l].resize(nmax_ + 1);
+      Psiz1[l].resize(nmax_ + 1);
+      Zetaz[l].resize(nmax_ + 1);
+      Zetaz1[l].resize(nmax_ + 1);
+    }
 
-  //*********************************************************************//
-  // Finally, we calculate the scattering coefficients (an and bn) and   //
-  // the angular functions (Pi and Tau). Note that for these arrays the  //
-  // first layer is 0 (zero), in future versions all arrays will follow  //
-  // this convention to save memory. (13 Nov, 2014)                      //
-  //*********************************************************************//
-  for (n = 0; n < nmax; n++) {
-    //********************************************************************//
-    //Expressions for calculating an and bn coefficients are not valid if //
-    //there is only one PEC layer (ie, for a simple PEC sphere).          //
-    //********************************************************************//
-    if (pl < (L - 1)) {
-      an[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
-      bn[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
-    } else {
-      an[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
-      bn[n] = PsiXL[n + 1]/ZetaXL[n + 1];
+    for (int l = 0; l < L; ++l) {
+      calcD1D3(z[l], D1z[l], D3z[l]);
+      calcD1D3(z1[l], D1z1[l], D3z1[l]);
+      calcPsiZeta(z[l], D1z[l], D3z[l], Psiz[l], Zetaz[l]);
+      calcPsiZeta(z1[l], D1z1[l], D3z1[l], Psiz1[l], Zetaz1[l]);
+    }
+    auto& m = refr_index_;
+    std::vector< std::complex<double> > m1(L);
+
+    for (int l = 0; l < L - 1; ++l) m1[l] = m[l + 1];
+    m1[L - 1] = std::complex<double> (1.0, 0.0);
+
+    // for (auto zz : m) printf ("m[i]=%g \n\n ", zz.real());
+    for (int l = L - 1; l >= 0; l--) {
+      for (int n = nmax_ - 2; n >= 0; n--) {
+        auto denomZeta = m1[l]*Zetaz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+        auto denomPsi = m1[l]*Psiz[l][n + 1]*(D1z[l][n + 1] - D3z[l][n + 1]);
+
+        auto T1 = anl_[l + 1][n]*Zetaz1[l][n + 1] - dnl_[l + 1][n]*Psiz1[l][n + 1];
+        auto T2 = bnl_[l + 1][n]*Zetaz1[l][n + 1] - cnl_[l + 1][n]*Psiz1[l][n + 1];
+
+        auto T3 = -D1z1[l][n + 1]*dnl_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*anl_[l + 1][n]*Zetaz1[l][n + 1];
+        auto T4 = -D1z1[l][n + 1]*cnl_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*bnl_[l + 1][n]*Zetaz1[l][n + 1];
+
+        // anl
+        anl_[l][n] = (D1z[l][n + 1]*m1[l]*T1 - m[l]*T3)/denomZeta;
+        // bnl
+        bnl_[l][n] = (D1z[l][n + 1]*m[l]*T2 - m1[l]*T4)/denomZeta;
+        // cnl
+        cnl_[l][n] = (D3z[l][n + 1]*m[l]*T2 - m1[l]*T4)/denomPsi;
+        // dnl
+        dnl_[l][n] = (D3z[l][n + 1]*m1[l]*T1 - m[l]*T3)/denomPsi;
+      }  // end of all n
+    }  // end of all l
+
+    // Check the result and change  an__0 and bn__0 for exact zero
+    for (int n = 0; n < nmax_; ++n) {
+      if (std::abs(anl_[0][n]) < 1e-10) anl_[0][n] = 0.0;
+      else throw std::invalid_argument("Unstable calculation of a__0_n!");
+      if (std::abs(bnl_[0][n]) < 1e-10) bnl_[0][n] = 0.0;
+      else throw std::invalid_argument("Unstable calculation of b__0_n!");
     }
-  }
 
-  return nmax;
-}
+    // for (int l = 0; l < L; ++l) {
+    //   printf("l=%d --> ", l);
+    //   for (int n = 0; n < nmax_ + 1; ++n) {
+    //         if (n < 20) continue;
+    //         printf("n=%d --> D1zn=%g, D3zn=%g, D1zn=%g, D3zn=%g || ",
+    //                n,
+    //                D1z[l][n].real(), D3z[l][n].real(),
+    //                D1z1[l][n].real(), D3z1[l][n].real());
+    //   }
+    //   printf("\n\n");
+    // }
+    // for (int l = 0; l < L; ++l) {
+    //   printf("l=%d --> ", l);
+    //   for (int n = 0; n < nmax_ + 1; ++n) {
+    //         printf("n=%d --> D1zn=%g, D3zn=%g, D1zn=%g, D3zn=%g || ",
+    //                n,
+    //                D1z[l][n].real(), D3z[l][n].real(),
+    //                D1z1[l][n].real(), D3z1[l][n].real());
+    //   }
+    //   printf("\n\n");
+    // }
+    //for (int i = 0; i < L + 1; ++i) {
+    //  printf("Layer =%d ---> \n", i);
+    //  for (int n = 0; n < nmax_; ++n) {
+    //                if (n < 20) continue;
+    //        printf(" || n=%d --> a=%g,%g b=%g,%g c=%g,%g d=%g,%g\n",
+    //               n,
+    //               anl_[i][n].real(), anl_[i][n].imag(),
+    //               bnl_[i][n].real(), bnl_[i][n].imag(),
+    //               cnl_[i][n].real(), cnl_[i][n].imag(),
+    //               dnl_[i][n].real(), dnl_[i][n].imag());
+    //  }
+    //  printf("\n\n");
+    //}
+    areIntCoeffsCalc_ = true;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
 
-//**********************************************************************************//
-// This function calculates the actual scattering parameters and amplitudes         //
-//                                                                                  //
-// 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]     //
-//   nTheta: Number of scattering angles                                            //
-//   Theta: Array containing all the scattering angles where the scattering         //
-//          amplitudes will be calculated                                           //
-//   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:                                                               //
-//   Qext: Efficiency factor for extinction                                         //
-//   Qsca: Efficiency factor for scattering                                         //
-//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
-//   Qbk: Efficiency factor for backscattering                                      //
-//   Qpr: Efficiency factor for the radiation pressure                              //
-//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
-//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
-//   S1, S2: Complex scattering amplitudes                                          //
-//                                                                                  //
-// Return value:                                                                    //
-//   Number of multipolar expansion terms used for the calculations                 //
-//**********************************************************************************//
 
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, 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 i, n, t;
-  std::vector<std::complex<double> > an, bn;
-  std::complex<double> Qbktmp;
-
-  // Calculate scattering coefficients
-  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
-
-  std::vector<double> Pi, Tau;
-  Pi.resize(nmax);
-  Tau.resize(nmax);
-
-  double x2 = x[L - 1]*x[L - 1];
-
-  // Initialize the scattering parameters
-  *Qext = 0;
-  *Qsca = 0;
-  *Qabs = 0;
-  *Qbk = 0;
-  Qbktmp = std::complex<double>(0.0, 0.0);
-  *Qpr = 0;
-  *g = 0;
-  *Albedo = 0;
-
-  // Initialize the scattering amplitudes
-  for (t = 0; t < nTheta; t++) {
-    S1[t] = std::complex<double>(0.0, 0.0);
-    S2[t] = std::complex<double>(0.0, 0.0);
-  }
-
-  // By using downward recurrence we avoid loss of precision due to float rounding errors
-  // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
-  //      http://en.wikipedia.org/wiki/Loss_of_significance
-  for (i = nmax - 2; i >= 0; i--) {
-    n = i + 1;
-    // Equation (27)
-    *Qext += (n + n + 1)*(an[i].real() + bn[i].real());
-    // Equation (28)
-    *Qsca += (n + n + 1)*(an[i].real()*an[i].real() + an[i].imag()*an[i].imag() + bn[i].real()*bn[i].real() + bn[i].imag()*bn[i].imag());
-    // Equation (29) TODO We must check carefully this equation. If we
-    // remove the typecast to double then the result changes. Which is
-    // the correct one??? Ovidio (2014/12/10) With cast ratio will
-    // give double, without cast (n + n + 1)/(n*(n + 1)) will be
-    // rounded to integer. Tig (2015/02/24)
-    *Qpr += ((n*(n + 2)/(n + 1))*((an[i]*std::conj(an[n]) + bn[i]*std::conj(bn[n])).real()) + ((double)(n + n + 1)/(n*(n + 1)))*(an[i]*std::conj(bn[i])).real());
-    // Equation (33)
-    Qbktmp = Qbktmp + (double)(n + n + 1)*(1 - 2*(n % 2))*(an[i]- bn[i]);
+  // ********************************************************************** //
+  // external scattering field = incident + scattered                       //
+  // BH p.92 (4.37), 94 (4.45), 95 (4.50)                                   //
+  // assume: medium is non-absorbing; refim = 0; Uabs = 0                   //
+  // ********************************************************************** //
 
-    //****************************************************//
-    // Calculate the scattering amplitudes (S1 and S2)    //
-    // Equations (25a) - (25b)                            //
-    //****************************************************//
-    for (t = 0; t < nTheta; t++) {
-      calcPiTau(nmax, Theta[t], Pi, Tau);
+  void MultiLayerMie::fieldExt(const double Rho, const double Phi, const double Theta,
+                               std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 
-      S1[t] += calc_S1(n, an[i], bn[i], Pi[i], Tau[i]);
-      S2[t] += calc_S2(n, an[i], bn[i], Pi[i], Tau[i]);
-    }
-  }
+    std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
+    std::vector<std::complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
+    std::vector<std::complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
+    std::vector<std::complex<double> > Ei(3, c_zero), Hi(3, c_zero), Es(3, c_zero), Hs(3, c_zero);
+    std::vector<std::complex<double> > jn(nmax_ + 1), jnp(nmax_ + 1), h1n(nmax_ + 1), h1np(nmax_ + 1);
+    std::vector<double> Pi(nmax_), Tau(nmax_);
 
-  *Qext = 2*(*Qext)/x2;                                 // Equation (27)
-  *Qsca = 2*(*Qsca)/x2;                                 // Equation (28)
-  *Qpr = *Qext - 4*(*Qpr)/x2;                           // Equation (29)
+    // Calculate spherical Bessel and Hankel functions
+    sbesjh(Rho, jn, jnp, h1n, h1np);
 
-  *Qabs = *Qext - *Qsca;                                // Equation (30)
-  *Albedo = *Qsca / *Qext;                              // Equation (31)
-  *g = (*Qext - *Qpr) / *Qsca;                          // Equation (32)
+    // Calculate angular functions Pi and Tau
+    calcPiTau(std::cos(Theta), Pi, Tau);
 
-  *Qbk = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2;    // Equation (33)
+    for (int n = 0; n < nmax_; n++) {
+      int n1 = n + 1;
+      double rn = static_cast<double>(n1);
 
-  return nmax;
-}
+      // using BH 4.12 and 4.50
+      calcSpherHarm(Rho, Phi, Theta, h1n[n1], h1np[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
 
-//**********************************************************************************//
-// This function is just a wrapper to call the full 'nMie' function with fewer      //
-// parameters, it is here mainly for compatibility with older versions of the       //
-// program. Also, you can use it if you neither have a PEC layer nor want to define //
-// any limit for the maximum number of terms.                                       //
-//                                                                                  //
-// Input parameters:                                                                //
-//   L: Number of layers                                                            //
-//   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]     //
-//   nTheta: Number of scattering angles                                            //
-//   Theta: Array containing all the scattering angles where the scattering         //
-//          amplitudes will be calculated                                           //
-//                                                                                  //
-// Output parameters:                                                               //
-//   Qext: Efficiency factor for extinction                                         //
-//   Qsca: Efficiency factor for scattering                                         //
-//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
-//   Qbk: Efficiency factor for backscattering                                      //
-//   Qpr: Efficiency factor for the radiation pressure                              //
-//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
-//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
-//   S1, S2: Complex scattering amplitudes                                          //
-//                                                                                  //
-// Return value:                                                                    //
-//   Number of multipolar expansion terms used for the calculations                 //
-//**********************************************************************************//
+      // scattered field: BH p.94 (4.45)
+      std::complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
+      for (int i = 0; i < 3; i++) {
+        Es[i] = Es[i] + En*(c_i*an_[n]*N3e1n[i] - bn_[n]*M3o1n[i]);
+        Hs[i] = Hs[i] + En*(c_i*bn_[n]*N3o1n[i] + an_[n]*M3e1n[i]);
+      }
+    }
 
-int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
-         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) {
+    // incident E field: BH p.89 (4.21); cf. p.92 (4.37), p.93 (4.38)
+    // basis unit vectors = er, etheta, ephi
+    std::complex<double> eifac = std::exp(std::complex<double>(0.0, Rho*std::cos(Theta)));
+    {
+      using std::sin;
+      using std::cos;
+      Ei[0] = eifac*sin(Theta)*cos(Phi);
+      Ei[1] = eifac*cos(Theta)*cos(Phi);
+      Ei[2] = -eifac*sin(Phi);
+    }
 
-  return nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
-}
+    // magnetic field
+    double hffact = 1.0/(cc_*mu_);
+    for (int i = 0; i < 3; i++) {
+      Hs[i] = hffact*Hs[i];
+    }
 
+    // incident H field: BH p.26 (2.43), p.89 (4.21)
+    std::complex<double> hffacta = hffact;
+    std::complex<double> hifac = eifac*hffacta;
+    {
+      using std::sin;
+      using std::cos;
+      Hi[0] = hifac*sin(Theta)*sin(Phi);
+      Hi[1] = hifac*cos(Theta)*sin(Phi);
+      Hi[2] = hifac*cos(Phi);
+    }
 
-//**********************************************************************************//
-// This function is just a wrapper to call the full 'nMie' function with fewer      //
-// parameters, it is useful if you want to include a PEC layer but not a limit      //
-// for the maximum number of terms.                                                 //
-//                                                                                  //
-// 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]     //
-//   nTheta: Number of scattering angles                                            //
-//   Theta: Array containing all the scattering angles where the scattering         //
-//          amplitudes will be calculated                                           //
-//                                                                                  //
-// Output parameters:                                                               //
-//   Qext: Efficiency factor for extinction                                         //
-//   Qsca: Efficiency factor for scattering                                         //
-//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
-//   Qbk: Efficiency factor for backscattering                                      //
-//   Qpr: Efficiency factor for the radiation pressure                              //
-//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
-//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
-//   S1, S2: Complex scattering amplitudes                                          //
-//                                                                                  //
-// Return value:                                                                    //
-//   Number of multipolar expansion terms used for the calculations                 //
-//**********************************************************************************//
+    for (int i = 0; i < 3; i++) {
+      // electric field E [V m - 1] = EF*E0
+      E[i] = Ei[i] + Es[i];
+      H[i] = Hi[i] + Hs[i];
+      // printf("ext E[%d]=%g",i,std::abs(E[i]));
+    }
+   }  // end of MultiLayerMie::fieldExt(...)
 
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         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) {
 
-  return nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
-}
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMie::fieldInt(const double Rho, const double Phi, const double Theta,
+                               std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 
-//**********************************************************************************//
-// This function is just a wrapper to call the full 'nMie' function with fewer      //
-// parameters, it is useful if you want to include a limit for the maximum number   //
-// of terms but not a PEC layer.                                                    //
-//                                                                                  //
-// Input parameters:                                                                //
-//   L: Number of layers                                                            //
-//   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]     //
-//   nTheta: Number of scattering angles                                            //
-//   Theta: Array containing all the scattering angles where the scattering         //
-//          amplitudes will be calculated                                           //
-//   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:                                                               //
-//   Qext: Efficiency factor for extinction                                         //
-//   Qsca: Efficiency factor for scattering                                         //
-//   Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)                    //
-//   Qbk: Efficiency factor for backscattering                                      //
-//   Qpr: Efficiency factor for the radiation pressure                              //
-//   g: Asymmetry factor (g = (Qext-Qpr)/Qsca)                                      //
-//   Albedo: Single scattering albedo (Albedo = Qsca/Qext)                          //
-//   S1, S2: Complex scattering amplitudes                                          //
-//                                                                                  //
-// Return value:                                                                    //
-//   Number of multipolar expansion terms used for the calculations                 //
-//**********************************************************************************//
+    std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
+    std::vector<std::complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
+    std::vector<std::complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
+    std::vector<std::complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
+    std::vector<std::complex<double> > El(3, c_zero), Ei(3, c_zero), Eic(3, c_zero), Hl(3, c_zero);
+    std::vector<std::complex<double> > jn(nmax_ + 1), jnp(nmax_ + 1), h1n(nmax_ + 1), h1np(nmax_ + 1);
+    std::vector<double> Pi(nmax_), Tau(nmax_);
 
-int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, 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 l = 0;  // Layer number
+    std::complex<double> ml;
 
-  return nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
-}
+    if (Rho > size_param_.back()) {
+      l = size_param_.size();
+      ml = c_one;
+    } else {
+      for (int i = 0; i < size_param_.size() - 1; ++i) {
+        if (size_param_[i] < Rho && Rho <= size_param_[i + 1]) {
+          l = i;
+        }
+      }
+      ml = refr_index_[l];
+    }
 
+//    for (int i = 0; i < size_param_.size(); i++) {
+//      printf("x[%i] = %g; m[%i] = %g, %g; ", i, size_param_[i], i, refr_index_[i].real(), refr_index_[i].imag());
+//    }
+//    printf("\nRho = %g; Phi = %g; Theta = %g; x[%i] = %g; m[%i] = %g, %g\n", Rho, Phi, Theta, l, size_param_[l], l, ml.real(), ml.imag());
 
-//**********************************************************************************//
-// This function calculates complex electric and magnetic field in the surroundings //
-// and inside (TODO) the particle.                                                  //
-//                                                                                  //
-// Input parameters:                                                                //
-//   L: Number of layers                                                            //
-//   pl: Index of PEC layer. If there is none just send 0 (zero)                    //
-//   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 0 (zero) and the function will calculate it.       //
-//   ncoord: Number of coordinate points                                            //
-//   Coords: Array containing all coordinates where the complex electric and        //
-//           magnetic fields will be calculated                                     //
-//                                                                                  //
-// Output parameters:                                                               //
-//   E, H: Complex electric and magnetic field at the provided coordinates          //
-//                                                                                  //
-// Return value:                                                                    //
-//   Number of multipolar expansion terms used for the calculations                 //
-//**********************************************************************************//
+    // Calculate spherical Bessel and Hankel functions
+    sbesjh(Rho*ml, jn, jnp, h1n, h1np);
 
-int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax, int ncoord, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
+    // Calculate angular functions Pi and Tau
+    calcPiTau(std::cos(Theta), Pi, Tau);
 
-  int i, c;
-  double Rho, Phi, Theta;
-  std::vector<std::complex<double> > an, bn;
+    for (int n = nmax_ - 1; n >= 0; n--) {
+      int n1 = n + 1;
+      double rn = static_cast<double>(n1);
 
-  // This array contains the fields in spherical coordinates
-  std::vector<std::complex<double> > Es, Hs;
-  Es.resize(3);
-  Hs.resize(3);
+      // using BH 4.12 and 4.50
+      calcSpherHarm(Rho, Phi, Theta, jn[n1], jnp[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
 
+      calcSpherHarm(Rho, Phi, Theta, h1n[n1], h1np[n], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
 
-  // Calculate scattering coefficients
-  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
+      // Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
+      std::complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
+      for (int i = 0; i < 3; i++) {
+        Ei[i] = Ei[i] + En*(M1o1n[i] - c_i*N1e1n[i]);
 
-  std::vector<double> Pi, Tau;
-  Pi.resize(nmax);
-  Tau.resize(nmax);
+        El[i] = El[i] + En*(cnl_[l][n]*M1o1n[i] - c_i*dnl_[l][n]*N1e1n[i]
+                      + c_i*anl_[l][n]*N3e1n[i] -     bnl_[l][n]*M3o1n[i]);
 
-  for (c = 0; c < ncoord; c++) {
-    // Convert to spherical coordinates
-    Rho = sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c] + Zp[c]*Zp[c]);
+        Hl[i] = Hl[i] + En*(-dnl_[l][n]*M1e1n[i] - c_i*cnl_[l][n]*N1o1n[i]
+                      +  c_i*bnl_[l][n]*N3o1n[i] +     anl_[l][n]*M3e1n[i]);
+      }
+    }  // end of for all n
+
+
+    // Debug field calculation outside the particle
+    if (l == size_param_.size()) {
+      // incident E field: BH p.89 (4.21); cf. p.92 (4.37), p.93 (4.38)
+      // basis unit vectors = er, etheta, ephi
+      std::complex<double> eifac = std::exp(std::complex<double>(0.0, Rho*std::cos(Theta)));
+      {
+        using std::sin;
+        using std::cos;
+        Eic[0] = eifac*sin(Theta)*cos(Phi);
+        Eic[1] = eifac*cos(Theta)*cos(Phi);
+        Eic[2] = -eifac*sin(Phi);
+      }
 
-    // Avoid convergence problems due to Rho too small
-    if (Rho < 1e-5) {
-      Rho = 1e-5;
+      printf("Rho = %g; Phi = %g; Theta = %g\n", Rho, Phi, Theta);
+      for (int i = 0; i < 3; i++) {
+        printf("Ei[%i] = %g, %g; Eic[%i] = %g, %g\n", i, Ei[i].real(), Ei[i].imag(), i, Eic[i].real(), Eic[i].imag());
+      }
     }
 
-    //If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
-    if (Rho == 0.0) {
-      Theta = 0.0;
-    } else {
-      Theta = acos(Zp[c]/Rho);
-    }
 
-    //If Xp=Yp=0 then Phi is undefined. Just set it to zero to zero to avoid problems
-    if ((Xp[c] == 0.0) and (Yp[c] == 0.0)) {
-      Phi = 0.0;
-    } else {
-      Phi = acos(Xp[c]/sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c]));
+    // magnetic field
+    double hffact = 1.0/(cc_*mu_);
+    for (int i = 0; i < 3; i++) {
+      Hl[i] = hffact*Hl[i];
     }
 
-    calcPiTau(nmax, Theta, Pi, Tau);
-
-    //*******************************************************//
-    // external scattering field = incident + scattered      //
-    // BH p.92 (4.37), 94 (4.45), 95 (4.50)                  //
-    // assume: medium is non-absorbing; refim = 0; Uabs = 0  //
-    //*******************************************************//
-
-    // Firstly the easiest case: the field outside the particle
-    if (Rho >= x[L - 1]) {
-      fieldExt(nmax, Rho, Phi, Theta, Pi, Tau, an, bn, Es, Hs);
-    } else {
-      // TODO, for now just set all the fields to zero
-      for (i = 0; i < 3; i++) {
-        Es[i] = std::complex<double>(0.0, 0.0);
-        Hs[i] = std::complex<double>(0.0, 0.0);
-      }
+    for (int i = 0; i < 3; i++) {
+      // electric field E [V m - 1] = EF*E0
+      E[i] = El[i];
+      H[i] = Hl[i];
     }
+   }  // end of MultiLayerMie::fieldInt(...)
+
+
+  //**********************************************************************************//
+  // This function calculates complex electric and magnetic field in the surroundings //
+  // and inside (TODO) the particle.                                                  //
+  //                                                                                  //
+  // Input parameters:                                                                //
+  //   L: Number of layers                                                            //
+  //   pl: Index of PEC layer. If there is none just send 0 (zero)                    //
+  //   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 0 (zero) and the function will calculate it.       //
+  //   ncoord: Number of coordinate points                                            //
+  //   Coords: Array containing all coordinates where the complex electric and        //
+  //           magnetic fields will be calculated                                     //
+  //                                                                                  //
+  // Output parameters:                                                               //
+  //   E, H: Complex electric and magnetic field at the provided coordinates          //
+  //                                                                                  //
+  // Return value:                                                                    //
+  //   Number of multipolar expansion terms used for the calculations                 //
+  //**********************************************************************************//
+  void MultiLayerMie::RunFieldCalculation() {
+    // Calculate external scattering coefficients an_ and bn_
+    ExtScattCoeffs();
+    // Calculate internal scattering coefficients anl_ and bnl_
+    IntScattCoeffs();
+
+//    for (int i = 0; i < nmax_; i++) {
+//      printf("a[%i] = %g, %g; b[%i] = %g, %g\n", i, an_[i].real(), an_[i].imag(), i, bn_[i].real(), bn_[i].imag());
+//    }
+
+    long total_points = coords_[0].size();
+    E_.resize(total_points);
+    H_.resize(total_points);
+    for (auto& f : E_) f.resize(3);
+    for (auto& f : H_) f.resize(3);
+
+    for (int point = 0; point < total_points; point++) {
+      const double& Xp = coords_[0][point];
+      const double& Yp = coords_[1][point];
+      const double& Zp = coords_[2][point];
+
+      // Convert to spherical coordinates
+      double Rho = std::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
+
+      // If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
+      double Theta = (Rho > 0.0) ? std::acos(Zp/Rho) : 0.0;
+
+      // If Xp=Yp=0 then Phi is undefined. Just set it to zero to avoid problems
+      double Phi = (Xp != 0.0 || Yp != 0.0) ? std::acos(Xp/std::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
+
+      // Avoid convergence problems due to Rho too small
+      if (Rho < 1e-5) Rho = 1e-5;
+
+      //*******************************************************//
+      // external scattering field = incident + scattered      //
+      // BH p.92 (4.37), 94 (4.45), 95 (4.50)                  //
+      // assume: medium is non-absorbing; refim = 0; Uabs = 0  //
+      //*******************************************************//
+
+      // This array contains the fields in spherical coordinates
+      std::vector<std::complex<double> > Es(3), Hs(3);
+
+      // Firstly the easiest case: the field outside the particle
+      if (Rho >= GetSizeParameter()) {
+        fieldInt(Rho, Phi, Theta, Es, Hs);
+//        fieldExt(Rho, Phi, Theta, Es, Hs);
+        //printf("\nFin E ext: %g,%g,%g   Rho=%g\n", std::abs(Es[0]), std::abs(Es[1]),std::abs(Es[2]), Rho);
+      } else {
+        fieldInt(Rho, Phi, Theta, Es, Hs);
+//        printf("\nFin E int: %g,%g,%g   Rho=%g\n", std::abs(Es[0]), std::abs(Es[1]),std::abs(Es[2]), Rho);
+      }
 
-    //Now, convert the fields back to cartesian coordinates
-    E[c][0] = std::sin(Theta)*std::cos(Phi)*Es[0] + std::cos(Theta)*std::cos(Phi)*Es[1] - std::sin(Phi)*Es[2];
-    E[c][1] = std::sin(Theta)*std::sin(Phi)*Es[0] + std::cos(Theta)*std::sin(Phi)*Es[1] + std::cos(Phi)*Es[2];
-    E[c][2] = std::cos(Theta)*Es[0] - std::sin(Theta)*Es[1];
+      //Now, convert the fields back to cartesian coordinates
+      {
+        using std::sin;
+        using std::cos;
+        E_[point][0] = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
+        E_[point][1] = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
+        E_[point][2] = cos(Theta)*Es[0] - sin(Theta)*Es[1];
+
+        H_[point][0] = sin(Theta)*cos(Phi)*Hs[0] + cos(Theta)*cos(Phi)*Hs[1] - sin(Phi)*Hs[2];
+        H_[point][1] = sin(Theta)*sin(Phi)*Hs[0] + cos(Theta)*sin(Phi)*Hs[1] + cos(Phi)*Hs[2];
+        H_[point][2] = cos(Theta)*Hs[0] - sin(Theta)*Hs[1];
+      }
+      //printf("Cart E: %g,%g,%g   Rho=%g\n", std::abs(Ex), std::abs(Ey),std::abs(Ez), Rho);
+    }  // end of for all field coordinates
 
-    H[c][0] = std::sin(Theta)*std::cos(Phi)*Hs[0] + std::cos(Theta)*std::cos(Phi)*Hs[1] - std::sin(Phi)*Hs[2];
-    H[c][1] = std::sin(Theta)*std::sin(Phi)*Hs[0] + std::cos(Theta)*std::sin(Phi)*Hs[1] + std::cos(Phi)*Hs[2];
-    H[c][2] = std::cos(Theta)*Hs[0] - std::sin(Theta)*Hs[1];
-  }
+  }  //  end of MultiLayerMie::RunFieldCalculation()
 
-  return nmax;
-}
+}  // end of namespace nmie

nmie.h

@@ -1,5 +1,6 @@
 //**********************************************************************************//
 //    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
 //                                                                                  //
 //    This file is part of scattnlay                                                //
 //                                                                                  //
@@ -25,34 +26,158 @@
 //**********************************************************************************//
 
 #define VERSION "0.3.1"
+#include <array>
 #include <complex>
+#include <cstdlib>
+#include <iostream>
 #include <vector>
 
-int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
-		        std::vector<std::complex<double> > &an, std::vector<std::complex<double> > &bn);
+namespace nmie {
+  int ScattCoeffs(const 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 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, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+  int nMie(const int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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 int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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 int L, 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, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+  int nField(const int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const int ncoord, const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
 
-int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
-         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);
+  class MultiLayerMie {
+   public:
+    // Run calculation
+    void RunMieCalculation();
+    void RunFieldCalculation();
 
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         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);
+    // Return calculation results
+    double GetQext();
+    double GetQsca();
+    double GetQabs();
+    double GetQbk();
+    double GetQpr();
+    double GetAsymmetryFactor();
+    double GetAlbedo();
+    std::vector<std::complex<double> > GetS1();
+    std::vector<std::complex<double> > GetS2();
 
-int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, 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);
+    std::vector<std::complex<double> > GetAn(){return an_;};
+    std::vector<std::complex<double> > GetBn(){return bn_;};
 
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, 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);
+    // Problem definition
+    // Add new layer
+    void AddNewLayer(double layer_size, std::complex<double> layer_index);
+    // Modify width of the layer
+    void SetLayerSize(std::vector<double> layer_size, int layer_position = 0);
+    // Modify refractive index of the layer
+    void SetLayerIndex(std::vector< std::complex<double> > layer_index, int layer_position = 0);
+    // Modify size of all layers
+    void SetLayersSize(const std::vector<double>& layer_size);
+    // Modify refractive index of all layers
+    void SetLayersIndex(const std::vector< std::complex<double> >& index);
+    // Modify scattering (theta) angles
+    void SetAngles(const std::vector<double>& angles);
+    // Modify coordinates for field calculation
+    void SetFieldCoords(const std::vector< std::vector<double> >& coords);
+    // Modify PEC layer
+    void SetPECLayer(int layer_position = 0);
 
-int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
-           int ncoord, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
-		   std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
+    // Set a fixed value for the maximun number of terms
+    void SetMaxTerms(int nmax);
+    // Get maximun number of terms
+    int GetMaxTerms() {return nmax_;};
 
+    // Clear layer information
+    void ClearLayers();
 
+    // Applied units requests
+    double GetSizeParameter();
+    double GetLayerWidth(int layer_position = 0);
+    std::vector<double> GetLayersSize();
+    std::vector<std::complex<double> > GetLayersIndex();
+    std::vector<std::array<double, 3> > GetFieldCoords();
+
+    std::vector<std::vector< std::complex<double> > > GetFieldE(){return E_;};   // {X[], Y[], Z[]}
+    std::vector<std::vector< std::complex<double> > > GetFieldH(){return H_;};
+  private:
+    void calcNstop();
+    void calcNmax(int first_layer);
+    void sbesjh(std::complex<double> z,
+                std::vector<std::complex<double> >& jn, std::vector<std::complex<double> >& jnp, 
+                std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np);
+    std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
+                                 std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                                 std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1);
+    std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
+                                 std::complex<double> PsiXL, std::complex<double> ZetaXL,
+                                 std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1);
+    std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
+                                 double Pi, double Tau);
+    std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,
+                                 double Pi, double Tau);
+    void calcPsiZeta(std::complex<double> x,
+                     std::vector<std::complex<double> > D1,
+                     std::vector<std::complex<double> > D3,
+                     std::vector<std::complex<double> >& Psi,
+                     std::vector<std::complex<double> >& Zeta);
+    void calcD1D3(std::complex<double> z,
+                  std::vector<std::complex<double> >& D1,
+                  std::vector<std::complex<double> >& D3);
+    void calcPiTau(const double& costheta,
+                   std::vector<double>& Pi, std::vector<double>& Tau);
+    void calcSpherHarm(const double Rho, const double Phi, const double Theta,
+                       const std::complex<double>& zn, const std::complex<double>& dzn,
+                       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 ExtScattCoeffs();
+    void IntScattCoeffs();
+
+    void fieldExt(const double Rho, const double Phi, const double Theta,
+                  std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
+
+    void fieldInt(const double Rho, const double Phi, const double Theta,
+                  std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
+
+    bool areIntCoeffsCalc_ = false;
+    bool areExtCoeffsCalc_ = false;
+    bool isMieCalculated_ = false;
+
+    // Size parameter for all layers
+    std::vector<double> size_param_;
+    // Refractive index for all layers
+    std::vector< std::complex<double> > refr_index_;
+    // Scattering angles for scattering pattern in radians
+    std::vector<double> theta_;
+    // Should be -1 if there is no PEC.
+    int PEC_layer_position_ = -1;
+
+    // with calcNmax(int first_layer);
+    int nmax_ = -1;
+    int nmax_preset_ = -1;
+    // Scattering coefficients
+    std::vector<std::complex<double> > an_, bn_;
+    std::vector< std::vector<double> > coords_;
+    // TODO: check if l index is reversed will lead to performance
+    // boost, if $a^(L+1)_n$ stored in anl_[n][0], $a^(L)_n$ in
+    // anl_[n][1] and so on...
+    // at the moment order is forward!
+    std::vector< std::vector<std::complex<double> > > anl_, bnl_, cnl_, dnl_;
+    /// Store result
+    double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
+    std::vector<std::vector< std::complex<double> > > E_, H_;  // {X[], Y[], Z[]}
+    // Mie efficinecy from each multipole channel.
+    std::vector<double> Qsca_ch_, Qext_ch_, Qabs_ch_, Qbk_ch_, Qpr_ch_;
+    std::vector<double> Qsca_ch_norm_, Qext_ch_norm_, Qabs_ch_norm_, Qbk_ch_norm_, Qpr_ch_norm_;
+    std::vector<std::complex<double> > S1_, S2_;
+
+    //Used constants
+    const double PI_=3.14159265358979323846;
+    // light speed [m s-1]
+    double const cc_ = 2.99792458e8;
+    // assume non-magnetic (MU=MU0=const) [N A-2]
+    double const mu_ = 4.0*PI_*1.0e-7;
+
+    //Temporary variables
+    std::vector<std::complex<double> > PsiZeta_;
+
+
+  };  // end of class MultiLayerMie
+
+}  // end of namespace nmie

py_nmie.cc

@@ -33,8 +33,8 @@
 // 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> >
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, int nmax,
+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,
 		 double S1r[], double S1i[], double S2r[], double S2i[]) {
 
@@ -43,7 +43,7 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
   S1.resize(nTheta);
   S2.resize(nTheta);
 
-  result = nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+  result = nmie::nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 
   for (i = 0; i < nTheta; i++) {
     S1r[i] = S1[i].real();
@@ -58,8 +58,8 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 // Same as nField 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 nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
-           int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
+int nField(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+           const int nCoords, std::vector<double>& Xp, std::vector<double>& Yp, std::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,7 +72,7 @@ int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double
     H[i].resize(3);
   }
 
-  result = nField(L, pl, x, m, nmax, nCoords, Xp, Yp, Zp, E, H);
+  result = nmie::nField(L, pl, x, m, nmax, nCoords, Xp, Yp, Zp, E, H);
 
   for (i = 0; i < nCoords; i++) {
     Erx[i] = E[i][0].real();

py_nmie.h

@@ -27,13 +27,13 @@
 #include <complex>
 #include <vector>
 
-int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, int nmax,
+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,
 		 double S1r[], double S1i[], double S2r[], double S2i[]);
 
-int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
-           int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
+int nField(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+           const int nCoords, std::vector<double>& Xp, std::vector<double>& Yp, std::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[]);
 

scattnlay.pyx

@@ -21,28 +21,16 @@
 #
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
 # distutils: language = c++
-# distutils: sources = nmie-wrapper.cc
+# distutils: sources = nmie.cc
+
 from __future__ import division
 import numpy as np
 cimport numpy as np
 from libcpp.vector cimport vector
 from libcpp.vector cimport complex
 
-# cdef extern from "<vector>" namespace "std":
-#     cdef cppclass vector[T]:
-#         cppclass iterator:
-#             T operator*()
-#             iterator operator++()
-#             bint operator==(iterator)
-#             bint operator!=(iterator)
-#         vector()
-#         void push_back(T&)
-#         T& operator[](int)
-#         T& at(int)
-#         iterator begin()
-#         iterator end()
-
 cdef inline double *npy2c(np.ndarray a):
     assert a.dtype == np.float64
 
@@ -51,112 +39,12 @@ cdef inline double *npy2c(np.ndarray a):
 
     # Return data pointer
     return <double *>(a.data)
-##############################################################################
-##############################################################################
-##############################################################################
-##############################################################################
-
-# cdef extern from "py_nmie.h":
-#     cdef int nMie(int L, int pl, vector[double] x, vector[double 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[double 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[])
-##############################################################################
-##############################################################################
-##############################################################################
-##############################################################################
-
-cdef extern from "nmie-wrapper.h"  namespace "nmie":
-    cdef int nMie_wrapper(int L, const vector[double] x, const vector[double complex] m , int nTheta, const vector[double] Theta, double *qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector[double complex] S1, vector[double complex] S2);
-
-
-##############################################################################
-##############################################################################
-##############################################################################
-##############################################################################
-# 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):
-#     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.float64_t, ndim = 1] Qext = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] Qabs = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] Qsca = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] Qbk = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] Qpr = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] g = np.zeros(x.shape[0], dtype = np.float64)
-#     cdef np.ndarray[np.float64_t, ndim = 1] Albedo = np.zeros(x.shape[0], dtype = np.float64)
-
-#     cdef np.ndarray[np.complex128_t, ndim = 2] S1 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128)
-#     cdef np.ndarray[np.complex128_t, ndim = 2] S2 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128)
-
-#     cdef np.ndarray[np.float64_t, ndim = 1] S1r
-#     cdef np.ndarray[np.float64_t, ndim = 1] S1i
-#     cdef np.ndarray[np.float64_t, ndim = 1] S2r
-#     cdef np.ndarray[np.float64_t, ndim = 1] S2i
-
-#     for i in range(x.shape[0]):
-#         S1r = np.zeros(theta.shape[0], dtype = np.float64)
-#         S1i = np.zeros(theta.shape[0], dtype = np.float64)
-#         S2r = np.zeros(theta.shape[0], dtype = np.float64)
-#         S2i = np.zeros(theta.shape[0], dtype = np.float64)
-
-#         terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i))
-
-#         S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C')
-#         S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C')
-
-#     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):
-#     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 = 3] E = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
-#     cdef np.ndarray[np.complex128_t, ndim = 3] H = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
-
-#     cdef np.ndarray[np.float64_t, ndim = 1] Erx
-#     cdef np.ndarray[np.float64_t, ndim = 1] Ery
-#     cdef np.ndarray[np.float64_t, ndim = 1] Erz
-#     cdef np.ndarray[np.float64_t, ndim = 1] Eix
-#     cdef np.ndarray[np.float64_t, ndim = 1] Eiy
-#     cdef np.ndarray[np.float64_t, ndim = 1] Eiz
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hrx
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hry
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hrz
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hix
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hiy
-#     cdef np.ndarray[np.float64_t, ndim = 1] Hiz
-
-#     for i in range(x.shape[0]):
-#         Erx = np.zeros(coords.shape[0], dtype = np.float64)
-#         Ery = np.zeros(coords.shape[0], dtype = np.float64)
-#         Erz = np.zeros(coords.shape[0], dtype = np.float64)
-#         Eix = np.zeros(coords.shape[0], dtype = np.float64)
-#         Eiy = np.zeros(coords.shape[0], dtype = np.float64)
-#         Eiz = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hrx = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hry = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hrz = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hix = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hiy = np.zeros(coords.shape[0], dtype = np.float64)
-#         Hiz = np.zeros(coords.shape[0], dtype = np.float64)
-
-#         terms[i] = nField(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, coords.shape[0], coords[:, 0].copy('C'), coords[:, 1].copy('C'), coords[:, 2].copy('C'), npy2c(Erx), npy2c(Ery), npy2c(Erz), npy2c(Eix), npy2c(Eiy), npy2c(Eiz), npy2c(Hrx), npy2c(Hry), npy2c(Hrz), npy2c(Hix), npy2c(Hiy), npy2c(Hiz))
-
-#         E[i] = np.vstack((Erx.copy('C') + 1.0j*Eix.copy('C'), Ery.copy('C') + 1.0j*Eiy.copy('C'), Erz.copy('C') + 1.0j*Eiz.copy('C'))).transpose()
-#         H[i] = np.vstack((Hrx.copy('C') + 1.0j*Hix.copy('C'), Hry.copy('C') + 1.0j*Hiy.copy('C'), Hrz.copy('C') + 1.0j*Hiz.copy('C'))).transpose()
-
-#     return terms, E, H
-##############################################################################
-##############################################################################
-##############################################################################
-##############################################################################
-def scattnlay_wrapper(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):
+
+cdef extern from "py_nmie.h":
+    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):
     cdef Py_ssize_t i
 
     cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int)
@@ -183,15 +71,53 @@ def scattnlay_wrapper(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.comple
         S2r = np.zeros(theta.shape[0], dtype = np.float64)
         S2i = np.zeros(theta.shape[0], dtype = np.float64)
 
-        terms[i] = nMie_wrapper(x.shape[1], x[i].copy('C'),
-                                m[i].copy('C'), theta.shape[0],
-                                theta.copy('C'), &Qext[i], &Qsca[i],
-                                &Qabs[i], &Qbk[i], &Qpr[i], &g[i],
-                                &Albedo[i],
-                                S1r, S2r)
+        terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i))
 
-        S1[i] = S1r.copy('C')
-        S2[i] = S2r.copy('C')
+        S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C')
+        S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C')
 
     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):
+    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 = 3] E = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
+    cdef np.ndarray[np.complex128_t, ndim = 3] H = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128)
+
+    cdef np.ndarray[np.float64_t, ndim = 1] Erx
+    cdef np.ndarray[np.float64_t, ndim = 1] Ery
+    cdef np.ndarray[np.float64_t, ndim = 1] Erz
+    cdef np.ndarray[np.float64_t, ndim = 1] Eix
+    cdef np.ndarray[np.float64_t, ndim = 1] Eiy
+    cdef np.ndarray[np.float64_t, ndim = 1] Eiz
+    cdef np.ndarray[np.float64_t, ndim = 1] Hrx
+    cdef np.ndarray[np.float64_t, ndim = 1] Hry
+    cdef np.ndarray[np.float64_t, ndim = 1] Hrz
+    cdef np.ndarray[np.float64_t, ndim = 1] Hix
+    cdef np.ndarray[np.float64_t, ndim = 1] Hiy
+    cdef np.ndarray[np.float64_t, ndim = 1] Hiz
+
+    for i in range(x.shape[0]):
+        Erx = np.zeros(coords.shape[0], dtype = np.float64)
+        Ery = np.zeros(coords.shape[0], dtype = np.float64)
+        Erz = np.zeros(coords.shape[0], dtype = np.float64)
+        Eix = np.zeros(coords.shape[0], dtype = np.float64)
+        Eiy = np.zeros(coords.shape[0], dtype = np.float64)
+        Eiz = np.zeros(coords.shape[0], dtype = np.float64)
+        Hrx = np.zeros(coords.shape[0], dtype = np.float64)
+        Hry = np.zeros(coords.shape[0], dtype = np.float64)
+        Hrz = np.zeros(coords.shape[0], dtype = np.float64)
+        Hix = np.zeros(coords.shape[0], dtype = np.float64)
+        Hiy = np.zeros(coords.shape[0], dtype = np.float64)
+        Hiz = np.zeros(coords.shape[0], dtype = np.float64)
+
+        terms[i] = nField(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, coords.shape[0], coords[:, 0].copy('C'), coords[:, 1].copy('C'), coords[:, 2].copy('C'), npy2c(Erx), npy2c(Ery), npy2c(Erz), npy2c(Eix), npy2c(Eiy), npy2c(Eiz), npy2c(Hrx), npy2c(Hry), npy2c(Hrz), npy2c(Hix), npy2c(Hiy), npy2c(Hiz))
+
+        E[i] = np.vstack((Erx.copy('C') + 1.0j*Eix.copy('C'), Ery.copy('C') + 1.0j*Eiy.copy('C'), Erz.copy('C') + 1.0j*Eiz.copy('C'))).transpose()
+        H[i] = np.vstack((Hrx.copy('C') + 1.0j*Hix.copy('C'), Hry.copy('C') + 1.0j*Hiy.copy('C'), Hrz.copy('C') + 1.0j*Hiz.copy('C'))).transpose()
+
+    return terms, E, H
+

setup.py

@@ -36,31 +36,26 @@ from distutils.core import setup
 from distutils.extension import Extension
 import numpy as np
 
-# setup(name = __mod__,
-#       version = __version__,
-#       description = __title__,
-#       long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
-# radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
-# resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
-# considerable improvement over previous implementations of similar algorithms. For details see: \
-# O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
-#       author = __author__,
-#       author_email = __email__,
-#       maintainer = __author__,
-#       maintainer_email = __email__,
-#       keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
-#       url = __url__,
-#       license = 'GPL',
-#       platforms = 'any',
-#       ext_modules = [Extension("scattnlay", ["nmie.cc", "nmie-wrapper.cc", "py_nmie.cc", "scattnlay.cc"], language = "c++", include_dirs = [np.get_include()])], 
-#       extra_compile_args=['-std=c++11']
-# )
+setup(name = __mod__,
+      version = __version__,
+      description = __title__,
+      long_description="""The Python version of scattnlay, a computer implementation of the algorithm for the calculation of electromagnetic \
+radiation scattering by a multilayered sphere developed by Yang. It has been shown that the program is effective, \
+resulting in very accurate values of scattering efficiencies for a wide range of size parameters, which is a \
+considerable improvement over previous implementations of similar algorithms. For details see: \
+O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
+      author = __author__,
+      author_email = __email__,
+      maintainer = __author__,
+      maintainer_email = __email__,
+      keywords = ['Mie scattering', 'Multilayered sphere', 'Efficiency factors', 'Cross-sections'],
+      url = __url__,
+      license = 'GPL',
+      platforms = 'any',
+      ext_modules = [Extension("scattnlay",
+                               ["nmie.cc", "py_nmie.cc", "scattnlay.cpp"],
+                               language = "c++",
+                               include_dirs = [np.get_include()])], 
+      extra_compile_args=['-std=c++11']
+)
 
-from distutils.core import setup
-from Cython.Build import cythonize
-setup(ext_modules = cythonize(
-           "scattnlay.pyx",                 # our Cython source
-           sources=["nmie-wrapper.cc"],  # additional source file(s)
-           language="c++",             # generate C++ code
-           extra_compile_args=['-std=c++11'],
-      ))

setup_cython.py

@@ -34,7 +34,7 @@ __url__ = 'http://scattering.sourceforge.net/'
 
 from distutils.core import setup
 from distutils.extension import Extension
-from Cython.Distutils import build_ext
+from Cython.Build import cythonize
 import numpy as np
 
 setup(name = __mod__,
@@ -53,7 +53,10 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
       url = __url__,
       license = 'GPL',
       platforms = 'any',
-      cmdclass = {'build_ext': build_ext},
-      ext_modules = [Extension("scattnlay", ["nmie.cc", "py_nmie.cc", "scattnlay.pyx"], language = "c++", include_dirs = [np.get_include()])]
+      ext_modules = cythonize("scattnlay.pyx",                                        # our Cython source
+                              sources = ["nmie.cc", "py_nmie.cc", "scattnlay.cpp"],   # additional source file(s)
+                              language = "c++",                                       # generate C++ code
+                              extra_compile_args = ['-std=c++11'],
+      )
 )
 

standalone.cc

@@ -36,7 +36,6 @@
 #include <time.h>
 #include <string.h>
 #include "nmie.h"
-#include "nmie-wrapper.h"
 
 const double PI=3.14159265358979323846;
 
@@ -63,8 +62,8 @@ int main(int argc, char *argv[]) {
     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 ...] "
-			  + "[-t ti tf nt] [-c comment]\n");
+                          + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] "
+                          + "[-t ti tf nt] [-c comment]\n");
     enum mode_states {read_L, read_x, read_mr, read_mi, read_ti, read_tf, read_nt, read_comment};
     // for (auto arg : args) std::cout<< arg <<std::endl;
     std::string comment;
@@ -92,43 +91,43 @@ int main(int argc, char *argv[]) {
 
       // Detecting new read mode (if it is a valid -key) 
       if (arg == "-l") {
-	mode = read_L;
-	continue;
+        mode = read_L;
+        continue;
       }
       if (arg == "-t") {
-	if ((mode != read_x) && (mode != read_comment))
-	  throw std::invalid_argument(std::string("Unfinished layer!\n")
-							 +error_msg);
-	mode = read_ti;
-	continue;
+        if ((mode != read_x) && (mode != read_comment))
+          throw std::invalid_argument(std::string("Unfinished layer!\n")
+                                                         +error_msg);
+        mode = read_ti;
+        continue;
       }
       if (arg == "-c") {
-	if ((mode != read_x) && (mode != read_nt))
-	  throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg);
-	mode = read_comment;
-	continue;
+        if ((mode != read_x) && (mode != read_nt))
+          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;
+        L = std::stoi(arg);
+        mode = read_x;
+        continue;
       }
       if (mode == read_x) {
-	x.push_back(std::stod(arg));
-	mode = read_mr;
-	continue;
+        x.push_back(std::stod(arg));
+        mode = read_mr;
+        continue;
       }
       if (mode == read_mr) {
-	tmp_mr = std::stod(arg);
-	mode = read_mi;
-	continue;
+        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;
+        m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
+        mode = read_x;
+        continue;
       }
       // if (strcmp(argv[i], "-l") == 0) {
       //   i++;
@@ -136,7 +135,7 @@ int main(int argc, char *argv[]) {
       //   x.resize(L);
       //   m.resize(L);
       //   if (argc < 3*(L + 1)) {
-      // 	  throw std::invalid_argument(error_msg);
+      //           throw std::invalid_argument(error_msg);
       //   } else {
       //     for (l = 0; l < L; l++) {
       //       i++;
@@ -147,23 +146,23 @@ int main(int argc, char *argv[]) {
       //     }
       //   }
       if (mode == read_ti) {
-	ti = std::stod(arg);
-	mode = read_tf;
-	continue;
+        ti = std::stod(arg);
+        mode = read_tf;
+        continue;
       }
       if (mode == read_tf) {
-	tf = std::stod(arg);
-	mode = read_nt;
-	continue;
+        tf = std::stod(arg);
+        mode = read_nt;
+        continue;
       }
       if (mode == read_nt) {
-	nt = std::stoi(arg);
+        nt = std::stoi(arg);
         Theta.resize(nt);
         S1.resize(nt);
         S2.resize(nt);
         S1w.resize(nt);
         S2w.resize(nt);
-	continue;
+        continue;
       }
       //} else if (strcmp(argv[i], "-t") == 0) {
         // i++;
@@ -177,23 +176,23 @@ int main(int argc, char *argv[]) {
         // S1.resize(nt);
         // S2.resize(nt);
       if (mode ==  read_comment) {
-	comment = arg;
+        comment = arg;
         has_comment = 1;
-	continue;
+        continue;
       }
       // } else if (strcmp(argv[i], "-c") == 0) {
       //   i++;
-      // 	comment = args[i];
+      //         comment = args[i];
       //   //strcpy(comment, argv[i]);
       //   has_comment = 1;
       // } else { i++; }
     }
     if ( (x.size() != m.size()) || (L != x.size()) ) 
       throw std::invalid_argument(std::string("Broken structure!\n")
-							 +error_msg);
+                                                         +error_msg);
     if ( (0 == m.size()) || ( 0 == x.size()) ) 
       throw std::invalid_argument(std::string("Empty structure!\n")
-							 +error_msg);
+                                                         +error_msg);
     
     if (nt < 0) {
       printf("Error reading Theta.\n");
@@ -211,7 +210,7 @@ int main(int argc, char *argv[]) {
 
 
     //nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-    nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
+    nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
    
 
     if (has_comment) {

tests/python/field.py

@@ -42,15 +42,17 @@ print(scattnlay.__file__)
 from scattnlay import fieldnlay
 import numpy as np
 
-x = np.ones((1, 1), dtype = np.float64)
-x[0, 0] = 1.
+x = np.ones((1, 2), dtype = np.float64)
+x[0, 0] = 2.0*np.pi*0.05/1.064
+x[0, 1] = 2.0*np.pi*0.06/1.064
 
-m = np.ones((1, 1), dtype = np.complex128)
-m[0, 0] = (0.05 + 2.070j)/1.46
+m = np.ones((1, 2), dtype = np.complex128)
+m[0, 0] = 1.53413/1.3205
+m[0, 1] = (0.565838 + 7.23262j)/1.3205
 
-npts = 101
+npts = 501
 
-scan = np.linspace(-3.0*x[0, 0], 3.0*x[0, 0], npts)
+scan = np.linspace(-4.0*x[0, 1], 4.0*x[0, 1], npts)
 
 coordX, coordY = np.meshgrid(scan, scan)
 coordX.resize(npts*npts)
@@ -85,12 +87,12 @@ try:
     scale_y = np.linspace(min(coordY), max(coordY), npts)
 
     # Define scale ticks
-    min_tick = max(0.5, min(min_tick, np.amin(edata)))
+    min_tick = max(0.01, min(min_tick, np.amin(edata)))
     max_tick = max(max_tick, np.amax(edata))
     scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
 
     # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
-    cax = ax.imshow(edata, interpolation = 'bicubic', cmap = cm.afmhot,
+    cax = ax.imshow(edata, interpolation = 'nearest', cmap = cm.afmhot,
                     origin = 'lower', vmin = min_tick, vmax = max_tick,
                     extent = (min(scale_x), max(scale_x), min(scale_y), max(scale_y)),
                     norm = LogNorm())
@@ -104,9 +106,21 @@ try:
     plt.xlabel('X')
     plt.ylabel('Y')
 
+    # This part draws the nanoshell
+    from matplotlib import patches
+
+    s1 = patches.Arc((0, 0), 2.0*x[0, 0], 2.0*x[0, 0], angle=0.0, zorder=2,
+                      theta1=0.0, theta2=360.0, linewidth=1, color='#00fa9a')
+    ax.add_patch(s1)
+
+    s2 = patches.Arc((0, 0), 2.0*x[0, 1], 2.0*x[0, 1], angle=0.0, zorder=2,
+                      theta1=0.0, theta2=360.0, linewidth=1, color='#00fa9a')
+    ax.add_patch(s2)
+    # End of drawing
+
     plt.draw()
 
-    # plt.show()
+    plt.show()
 
     plt.clf()
     plt.close()

tests/python/lfield.py

@@ -37,11 +37,13 @@
 from scattnlay import fieldnlay
 import numpy as np
 
-x = np.ones((1, 1), dtype = np.float64)
-x[0, 0] = 1.
+x = np.ones((1, 2), dtype = np.float64)
+x[0, 0] = 2.0*np.pi*0.05/1.064
+x[0, 1] = 2.0*np.pi*0.06/1.064
 
-m = np.ones((1, 1), dtype = np.complex128)
-m[0, 0] = (0.05 + 2.070j)/1.46
+m = np.ones((1, 2), dtype = np.complex128)
+m[0, 0] = 1.53413/1.3205
+m[0, 1] = (0.565838 + 7.23262j)/1.3205
 
 nc = 1001
 
@@ -49,13 +51,19 @@ coordX = np.zeros((nc, 3), dtype = np.float64)
 coordY = np.zeros((nc, 3), dtype = np.float64)
 coordZ = np.zeros((nc, 3), dtype = np.float64)
 
-scan = np.linspace(-3.0*x[0, 0], 3.0*x[0, 0], nc)
+scan = np.linspace(-10.0*x[0, 1], 10.0*x[0, 1], nc)
 one = np.ones(nc, dtype = np.float64)
 
 coordX[:, 0] = scan
 coordY[:, 1] = scan
 coordZ[:, 2] = scan
 
+from scattnlay import scattnlay
+print "\nscattnlay"
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
+print "Results: ", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
+
+print "\nfieldnlay"
 terms, Ex, Hx = fieldnlay(x, m, coordX)
 terms, Ey, Hy = fieldnlay(x, m, coordY)
 terms, Ez, Hz = fieldnlay(x, m, coordZ)

tests/python/pfield.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#
+#    This file is part of python-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 along three 
+# points, for an spherical silver nanoparticle embedded in glass.
+
+# 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.
+
+from scattnlay import fieldnlay
+import numpy as np
+
+x = np.ones((1, 2), dtype = np.float64)
+x[0, 0] = 2.0*np.pi*0.05/1.064
+x[0, 1] = 2.0*np.pi*0.06/1.064
+
+m = np.ones((1, 2), dtype = np.complex128)
+m[0, 0] = 1.53413/1.3205
+m[0, 1] = (0.565838 + 7.23262j)/1.3205
+
+coord = np.zeros((3, 3), dtype = np.float64)
+coord[0, 0] = x[0, 0]/2.0
+coord[1, 0] = (x[0, 0] + x[0, 1])/2.0
+coord[2, 0] = 1.5*x[0, 1]
+
+terms, E, H = fieldnlay(x, m, coord)
+
+Er = np.absolute(E)
+
+# |E|/|Eo|
+Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+
+print x
+print m
+print np.vstack((coord[:, 0], Eh)).transpose()
+