The calculation of the electric field is almost done.

nmie.cc

@@ -87,104 +87,248 @@ int Nmax(int L, int fl, int pl,
 }
 
 //**********************************************************************************//
-<<<<<<< HEAD
-// This function calculates the spherical Bessel functions (jn and yn) for a given  //
-// real value r. See pag. 87 B&H.                                                   //
+// 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:                                                                //
-//   r: Real argument to evaluate jn and yn                                         //
-//   n_max: Maximum number of terms to calculate jn and yn                          //
+//   z: Real argument to evaluate jn and h1n                                        //
+//   nmax: Maximum number of terms to calculate jn and h1n                          //
 //                                                                                  //
 // Output parameters:                                                               //
-//   jn, yn: Spherical Bessel functions (double)                                    //
+//   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 sphericalBessel(double r, int n_max, std::vector<double> &j, std::vector<double> &y) {
+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, xdb, 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 (n_max >= 1) {
-    j[0] = sin(r)/r;
-    y[0] = -cos(r)/r;
+    if (absc < accur) {
+      // We have obtained the desired accuracy
+      break;
+    }
   }
 
-  if (n_max >= 2) {
-    j[1] = sin(r)/r/r - cos(r)/r;
-    y[1] = -cos(r)/r/r - sin(r)/r;
+  if (absc > accur) {
+    // We were not able to obtain the desired accuracy
+    return -2;
   }
 
-  for (n = 2; n < n_max; n++) {
-    j[n] = double(n + n + 1)*j[n - 1]/r - j[n - 2];
-    y[n] = double(n + n + 1)*y[n - 1]/r - h[n - 2];
-=======
-// This function calculates the spherical Bessel functions (jn and hn) for a given  //
-// value of z.                                                                      //
-//                                                                                  //
-// Input parameters:                                                                //
-//   z: Real argument to evaluate jn and hn                                         //
-//   n_max: Maximum number of terms to calculate jn and hn                          //
-//                                                                                  //
-// Output parameters:                                                               //
-//   jn, hn: Spherical Bessel functions (complex)                                   //
-//**********************************************************************************//
-void sphericalBessel(std::complex<double> r, int n_max, std::vector<std::complex<double> > &j, std::vector<std::complex<double> > &h) {
-  int n;
+  jn[nmax - 1] = tm30;
+  jnp[nmax - 1] = f*jn[nmax - 1];
 
-  j[0] = sin(r)/r;
-  j[1] = sin(r)/r/r - cos(r)/r;
-  h[0] = -cos(r)/r;
-  h[1] = -cos(r)/r/r - sin(r)/r;
-  for (n = 2; n < n_max; n++) {
-    j[n] = double(n + n + 1)*j[n - 1]/r - j[n - 2];
-    h[n] = double(n + n + 1)*h[n - 1]/r - h[n - 2];
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
+  // 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])) {
+        xdb = z;
+        h1nldb = h1n[n];
+        h1nbdb = h1n[n - 1];
+      }
+
+      pl += zi;
+
+      h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
+    }
+  }
+
+  // success
+  return 0;
 }
 
 //**********************************************************************************//
-<<<<<<< HEAD
-// This function calculates the spherical Hankel functions (h1n and h2n) for a      //
-// given real value r. See eqs. (4.13) and (4.14), pag. 87 B&H.                     //
+// 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:                                                                //
-//   r: Real argument to evaluate h1n and h2n                                       //
-//   n_max: Maximum number of terms to calculate h1n and h2n                        //
+//   z: Complex argument to evaluate bj, by and bd                                  //
+//   nmax: Maximum number of terms to calculate bj, by and bd                       //
 //                                                                                  //
 // Output parameters:                                                               //
-//   h1n, h2n: Spherical Hankel functions (complex)                                 //
+//   bj, by: Spherical Bessel functions                                             //
+//   bd: Logarithmic derivative                                                     //
 //**********************************************************************************//
-void sphericalHankel(double r, int n_max, std::vector<std::complex<double> > &h1, std::vector<std::complex<double> > &h2) {
-=======
-// This function calculates the spherical Bessel functions (jn and hn) for a given  //
-// value of r.                                                                      //
-//                                                                                  //
-// Input parameters:                                                                //
-//   r: Real argument to evaluate jn and hn                                         //
-//   n_max: Maximum number of terms to calculate jn and hn                          //
-//                                                                                  //
-// Output parameters:                                                               //
-//   jn, hn: Spherical Bessel functions (double)                                    //
-//**********************************************************************************//
-void sphericalBessel(double r, int n_max, std::vector<double> &j, std::vector<double> &h) {
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
-  int n;
-  std::complex<double> j, y;
-  j.resize(n_max);
-  h.resize(n_max);
-
-<<<<<<< HEAD
-  sphericalBessel(r, n_max, j, y);
-
-  for (n = 0; n < n_max; n++) {
-    h1[n] = std::complex<double> (j[n], y[n]);
-    h2[n] = std::complex<double> (j[n], -y[n]);
-=======
-  j[0] = sin(r)/r;
-  j[1] = sin(r)/r/r - cos(r)/r;
-  h[0] = -cos(r)/r;
-  h[1] = -cos(r)/r/r - sin(r)/r;
-  for (n = 2; n < n_max; n++) {
-    j[n] = double(n + n + 1)*j[n - 1]/r - j[n - 2];
-    h[n] = double(n + n + 1)*h[n - 1]/r - h[n - 2];
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
+void sphericalBessel(std::complex<double> z, int nmax, std::vector<double>& bj, std::vector<double>& by, std::vector<double>& bd) {
+
+    std::vector<std::complex<double> > jn, jnp, h1n, h1np;
+    jn.resize(nmax);
+    jnp.resize(nmax);
+    h1n.resize(nmax);
+    h1np.resize(nmax);
+
+    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
+int 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;
+  double rn = 0.0;
+  std::complex<double> 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);
+  by.resize(nmax);
+  bd.resize(nmax);
+
+  // Calculate spherical Bessel and Hankel functions
+  sphericalBessel(Rho, nmax, bj, by, bd);
+
+  for (n = 0; n < nmax; n++) {
+    rn = double(n + 1);
+
+    zn = bj[n] + std::complex<double>(0.0, 1.0)*by[n];
+    xxip = Rho*(bj[n] + std::complex<double>(0.0, 1.0)*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(std::complex<double>(0.0, 1.0), rn)*(2.0*rn + 1.0)/(rn*(rn + 1.0));
+    for (i = 0; i < 3; i++) {
+      Es[i] = Es[i] + encap*(std::complex<double>(0.0, 1.0)*an[n]*vn3e1n[i] - bn[n]*vm3o1n[i]);
+      Hs[i] = Hs[i] + encap*(std::complex<double>(0.0, 1.0)*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, 1.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
+  std::complex<double> hffact = std::complex<double>(refmed, 0.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];
   }
 }
 
@@ -232,23 +376,23 @@ std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double
 //                                                                                  //
 // Input parameters:                                                                //
 //   x: Real argument to evaluate Psi and Zeta                                      //
-//   n_max: Maximum number of terms to calculate 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 n_max,
+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) {
+		         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 <= n_max; n++) {
+  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]);
   }
@@ -261,29 +405,29 @@ void calcPsiZeta(double x, int n_max,
 //                                                                                  //
 // Input parameters:                                                                //
 //   z: Complex argument to evaluate D1 and D3                                      //
-//   n_max: Maximum number of terms to calculate 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 n_max,
-		      std::vector<std::complex<double> > &D1,
-		      std::vector<std::complex<double> > &D3) {
+void calcD1D3(std::complex<double> z, int nmax,
+		      std::vector<std::complex<double> >& D1,
+		      std::vector<std::complex<double> >& D3) {
 
   int n;
   std::vector<std::complex<double> > PsiZeta;
-  PsiZeta.resize(n_max + 1);
+  PsiZeta.resize(nmax + 1);
 
   // Downward recurrence for D1 - equations (16a) and (16b)
-  D1[n_max] = std::complex<double>(0.0, 0.0);
-  for (n = n_max; n > 0; n--) {
+  D1[nmax] = std::complex<double>(0.0, 0.0);
+  for (n = nmax; n > 0; n--) {
     D1[n - 1] = double(n)/z - 1.0/(D1[n] + double(n)/z);
   }
 
   // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
   PsiZeta[0] = 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 <= n_max; n++) {
+  for (n = 1; n <= nmax; n++) {
     PsiZeta[n] = PsiZeta[n - 1]*(double(n)/z - D1[n - 1])*(double(n)/z- D3[n - 1]);
     D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta[n];
   }
@@ -294,7 +438,7 @@ void calcD1D3(std::complex<double> z, int n_max,
 // Equations (26a) - (26c)                                                          //
 //                                                                                  //
 // Input parameters:                                                                //
-//   n_max: Maximum number of terms to calculate Pi and Tau                         //
+//   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                                           //
@@ -302,25 +446,25 @@ void calcD1D3(std::complex<double> z, int n_max,
 // Output parameters:                                                               //
 //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
 //**********************************************************************************//
-void calcPiTau(int n_max, int nTheta, std::vector<double> Theta,
-		       std::vector< std::vector<double> > &Pi,
-		       std::vector< std::vector<double> > &Tau) {
+void calcPiTau(int nmax, int nTheta, std::vector<double> Theta,
+		       std::vector< std::vector<double> >& Pi,
+		       std::vector< std::vector<double> >& Tau) {
 
   int n, t;
-  for (n = 0; n < n_max; n++) {
-    //****************************************************//
-    // Equations (26a) - (26c)                            //
-    //****************************************************//
-    for (t = 0; t < nTheta; t++) {
+  //****************************************************//
+  // Equations (26a) - (26c)                            //
+  //****************************************************//
+  for (t = 0; t < nTheta; t++) {
+    for (n = 0; n < nmax; n++) {
       if (n == 0) {
         // Initialize Pi and Tau
-        Pi[n][t] = 1.0;
-        Tau[n][t] = (n + 1)*cos(Theta[t]);
+        Pi[t][n] = 1.0;
+        Tau[t][n] = (n + 1)*cos(Theta[t]);
       } else {
         // Calculate the actual values
-        Pi[n][t] = ((n == 1) ? ((n + n + 1)*cos(Theta[t])*Pi[n - 1][t]/n)
-                             : (((n + n + 1)*cos(Theta[t])*Pi[n - 1][t] - (n + 1)*Pi[n - 2][t])/n));
-        Tau[n][t] = (n + 1)*cos(Theta[t])*Pi[n][t] - (n + 2)*Pi[n - 1][t];
+        Pi[t][n] = ((n == 1) ? ((n + n + 1)*cos(Theta[t])*Pi[t][n - 1]/n)
+                             : (((n + n + 1)*cos(Theta[t])*Pi[t][n - 1] - (n + 1)*Pi[t][n - 2])/n));
+        Tau[t][n] = (n + 1)*cos(Theta[t])*Pi[t][n] - (n + 2)*Pi[t][n - 1];
       }
     }
   }
@@ -335,7 +479,7 @@ void calcPiTau(int n_max, int nTheta, std::vector<double> Theta,
 //   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]     //
-//   n_max: 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 used if you know what you are doing, otherwise set   //
 //          this parameter to -1 and the function will calculate it.                //
 //                                                                                  //
@@ -345,8 +489,8 @@ void calcPiTau(int n_max, int nTheta, std::vector<double> Theta,
 // 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 n_max,
-		        std::vector<std::complex<double> > &an, std::vector<std::complex<double> > &bn) {
+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 //
@@ -355,8 +499,8 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
 
   int fl = (pl > 0) ? pl : 0;
 
-  if (n_max <= 0) {
-    n_max = Nmax(L, fl, pl, x, m);
+  if (nmax <= 0) {
+    nmax = Nmax(L, fl, pl, x, m);
   }
 
   std::complex<double> z1, z2;
@@ -391,35 +535,35 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
   Hb.resize(L);
 
   for (l = 0; l < L; l++) {
-    D1_mlxl[l].resize(n_max + 1);
-    D1_mlxlM1[l].resize(n_max + 1);
+    D1_mlxl[l].resize(nmax + 1);
+    D1_mlxlM1[l].resize(nmax + 1);
 
-    D3_mlxl[l].resize(n_max + 1);
-    D3_mlxlM1[l].resize(n_max + 1);
+    D3_mlxl[l].resize(nmax + 1);
+    D3_mlxlM1[l].resize(nmax + 1);
 
-    Q[l].resize(n_max + 1);
+    Q[l].resize(nmax + 1);
 
-    Ha[l].resize(n_max);
-    Hb[l].resize(n_max);
+    Ha[l].resize(nmax);
+    Hb[l].resize(nmax);
   }
 
-  an.resize(n_max);
-  bn.resize(n_max);
+  an.resize(nmax);
+  bn.resize(nmax);
 
   std::vector<std::complex<double> > D1XL, D3XL;
-  D1XL.resize(n_max + 1);
-  D3XL.resize(n_max + 1);
+  D1XL.resize(nmax + 1);
+  D3XL.resize(nmax + 1);
 
 
   std::vector<std::complex<double> > PsiXL, ZetaXL;
-  PsiXL.resize(n_max + 1);
-  ZetaXL.resize(n_max + 1);
+  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 <= n_max; n++) {
+    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);
     }
@@ -427,13 +571,13 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
     z1 = x[fl]* m[fl];
 
     // Calculate D1 and D3
-    calcD1D3(z1, n_max, D1_mlxl[fl], D3_mlxl[fl]);
+    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 < n_max; n++) {
+  for (n = 0; n < nmax; n++) {
     Ha[fl][n] = D1_mlxl[fl][n + 1];
     Hb[fl][n] = D1_mlxl[fl][n + 1];
   }
@@ -449,10 +593,10 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
     z2 = x[l - 1]*m[l];
 
     //Calculate D1 and D3 for z1
-    calcD1D3(z1, n_max, D1_mlxl[l], D3_mlxl[l]);
+    calcD1D3(z1, nmax, D1_mlxl[l], D3_mlxl[l]);
 
     //Calculate D1 and D3 for z2
-    calcD1D3(z2, n_max, D1_mlxlM1[l], D3_mlxlM1[l]);
+    calcD1D3(z2, nmax, D1_mlxlM1[l], D3_mlxlM1[l]);
 
     //*********************************************//
     //Calculate Q, Ha and Hb in the layers fl+1..L //
@@ -463,7 +607,7 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
     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 <= n_max; n++) {
+    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]);
 
@@ -471,7 +615,7 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
     }
 
     // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
-    for (n = 1; n <= n_max; n++) {
+    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];
@@ -511,10 +655,10 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
   //**************************************//
 
   // Calculate D1XL and D3XL
-  calcD1D3(x[L - 1], n_max, D1XL, D3XL);
+  calcD1D3(x[L - 1], nmax, D1XL, D3XL);
 
   // Calculate PsiXL and ZetaXL
-  calcPsiZeta(x[L - 1], n_max, D1XL, D3XL, PsiXL, ZetaXL);
+  calcPsiZeta(x[L - 1], nmax, D1XL, D3XL, PsiXL, ZetaXL);
 
   //*********************************************************************//
   // Finally, we calculate the scattering coefficients (an and bn) and   //
@@ -522,7 +666,7 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
   // first layer is 0 (zero), in future versions all arrays will follow  //
   // this convention to save memory. (13 Nov, 2014)                      //
   //*********************************************************************//
-  for (n = 0; n < n_max; n++) {
+  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).          //
@@ -536,7 +680,7 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
     }
   }
 
-  return n_max;
+  return nmax;
 }
 
 //**********************************************************************************//
@@ -550,7 +694,7 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
 //   nTheta: Number of scattering angles                                            //
 //   Theta: Array containing all the scattering angles where the scattering         //
 //          amplitudes will be calculated                                           //
-//   n_max: 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 used if you know what you are doing, otherwise set   //
 //          this parameter to -1 and the function will calculate it                 //
 //                                                                                  //
@@ -569,27 +713,26 @@ int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<d
 //**********************************************************************************//
 
 int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, int n_max,
+         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> >& 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
-  n_max = ScattCoeffs(L, pl, x, m, n_max, an, bn);
-
-  std::vector< std::vector<double> > Pi;
-  Pi.resize(n_max);
-  std::vector< std::vector<double> > Tau;
-  Tau.resize(n_max);
-  for (n = 0; n < n_max; n++) {
-    Pi[n].resize(nTheta);
-    Tau[n].resize(nTheta);
+  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
+
+  std::vector< std::vector<double> > Pi, Tau;
+  Pi.resize(nTheta);
+  Tau.resize(nTheta);
+  for (t = 0; t < nTheta; t++) {
+    Pi[t].resize(nmax);
+    Tau[t].resize(nmax);
   }
 
-  calcPiTau(n_max, nTheta, Theta, Pi, Tau);
+  calcPiTau(nmax, nTheta, Theta, Pi, Tau);
 
   double x2 = x[L - 1]*x[L - 1];
 
@@ -612,7 +755,7 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
   // 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 = n_max - 2; i >= 0; i--) {
+  for (i = nmax - 2; i >= 0; i--) {
     n = i + 1;
     // Equation (27)
     *Qext += (n + n + 1)*(an[i].real() + bn[i].real());
@@ -629,8 +772,8 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
     // Equations (25a) - (25b)                            //
     //****************************************************//
     for (t = 0; t < nTheta; t++) {
-      S1[t] += calc_S1(n, an[i], bn[i], Pi[i][t], Tau[i][t]);
-      S2[t] += calc_S2(n, an[i], bn[i], Pi[i][t], Tau[i][t]);
+      S1[t] += calc_S1(n, an[i], bn[i], Pi[t][i], Tau[t][i]);
+      S2[t] += calc_S2(n, an[i], bn[i], Pi[t][i], Tau[t][i]);
     }
   }
 
@@ -644,7 +787,7 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 
   *Qbk = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2;    // Equation (33)
 
-  return n_max;
+  return nmax;
 }
 
 //**********************************************************************************//
@@ -678,7 +821,7 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 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) {
+         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);
 }
@@ -715,7 +858,7 @@ int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
 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) {
+         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);
 }
@@ -732,7 +875,7 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 //   nTheta: Number of scattering angles                                            //
 //   Theta: Array containing all the scattering angles where the scattering         //
 //          amplitudes will be calculated                                           //
-//   n_max: 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 used if you know what you are doing, otherwise set   //
 //          this parameter to -1 and the function will calculate it                 //
 //                                                                                  //
@@ -751,11 +894,11 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
 //**********************************************************************************//
 
 int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
-         int nTheta, std::vector<double> Theta, int n_max,
+         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> >& S1, std::vector<std::complex<double> >& S2) {
 
-  return nMie(L, -1, x, m, nTheta, Theta, n_max, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
+  return nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 }
 
 
@@ -768,10 +911,10 @@ int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
 //   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]     //
-//   n_max: 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 used if you know what you are doing, otherwise set   //
 //          this parameter to 0 (zero) and the function will calculate it.          //
-//   nCoords: Number of coordinate points                                           //
+//   ncoord: Number of coordinate points                                           //
 //   Coords: Array containing all coordinates where the complex electric and        //
 //           magnetic fields will be calculated                                     //
 //                                                                                  //
@@ -782,44 +925,22 @@ int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
 //   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 n_max,
-           int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
-		   std::vector<std::complex<double> > &E, std::vector<std::complex<double> >  &H)  {
-<<<<<<< HEAD
-=======
-
-  int i, n, c;
-/*  double **Pi, **Tau;
-  complex *an, *bn;
-  double *Rho = (double *) malloc(nCoords*sizeof(double));
-  double *Phi = (double *) malloc(nCoords*sizeof(double));
-  double *Theta = (double *) malloc(nCoords*sizeof(double));
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
+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, n, c;
   std::vector<std::complex<double> > an, bn;
 
-<<<<<<< HEAD
   // Calculate scattering coefficients
-=======
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
-  n_max = ScattCoeffs(L, pl, x, m, n_max, an, bn);
-
-  std::vector< std::vector<double> > Pi, Tau;
-  Pi.resize(n_max);
-  Tau.resize(n_max);
-  for (n = 0; n < n_max; n++) {
-    Pi[n].resize(nCoords);
-    Tau[n].resize(nCoords);
-  }
+  nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
 
   std::vector<double> Rho, Phi, Theta;
-  Rho.resize(nCoords);
-  Phi.resize(nCoords);
-  Theta.resize(nCoords);
+  Rho.resize(ncoord);
+  Phi.resize(ncoord);
+  Theta.resize(ncoord);
 
-  for (c = 0; c < nCoords; c++) {
-<<<<<<< HEAD
+  for (c = 0; c < ncoord; c++) {
     // Convert to spherical coordinates
     Rho = sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c] + Zp[c]*Zp[c]);
     if (Rho < 1e-3) {
@@ -828,68 +949,33 @@ int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double
     Phi = acos(Xp[c]/sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c]));
     Theta = acos(Xp[c]/Rho[c]);
   }
-=======
-    E[c] = C_ZERO;
-    H[c] = C_ZERO;
-  }*/
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
-
-  calcPiTau(n_max, nCoords, Theta, Pi, Tau);
-
-<<<<<<< HEAD
-  std::vector<double > j, y;
-  std::vector<std::complex<double> > h1, h2;
-  j.resize(n_max);
-  y.resize(n_max);
-  h1.resize(n_max);
-  h2.resize(n_max);
-
-  for (c = 0; c < nCoords; c++) {
+
+  std::vector< std::vector<double> > Pi, Tau;
+  Pi.resize(ncoord);
+  Tau.resize(ncoord);
+  for (c = 0; c < ncoord; c++) {
+    Pi[c].resize(nmax);
+    Tau[c].resize(nmax);
+  }
+
+  calcPiTau(nmax, ncoord, Theta, Pi, Tau);
+
+  for (c = 0; c < ncoord; c++) {
     //*******************************************************//
     // 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 spherical Bessel and Hankel functions
-    sphericalBessel(Rho, n_max, j, y);
-    sphericalHankel(Rho, n_max, h1, h2);
-
     // Initialize the fields
     E[c] = std::complex<double>(0.0, 0.0);
     H[c] = std::complex<double>(0.0, 0.0);
 
     // Firstly the easiest case, we want the field outside the particle
     if (Rho >= x[L - 1]) {
-      
+      fieldExt(nmax, Rho, Phi, Theta, Pi[c], Tau[c], an, bn, E[c], H[c]);
     }
-=======
-  // Firstly the easiest case, we want the field outside the particle
-//  if (Rho[c] >= x[L - 1]) {
-//  }
-  for (i = 1; i < (n_max - 1); i++) {
-//    n = i - 1;
-/*    // Equation (27)
-    *Qext = *Qext + (double)(n + n + 1)*(an[i].r + bn[i].r);
-    // Equation (28)
-    *Qsca = *Qsca + (double)(n + n + 1)*(an[i].r*an[i].r + an[i].i*an[i].i + bn[i].r*bn[i].r + bn[i].i*bn[i].i);
-    // Equation (29)
-    *Qpr = *Qpr + ((n*(n + 2)/(n + 1))*((Cadd(Cmul(an[i], std::conj(an[n])), Cmul(bn[i], std::conj(bn[n])))).r) + ((double)(n + n + 1)/(n*(n + 1)))*(Cmul(an[i], std::conj(bn[i])).r));
-
-    // Equation (33)
-    Qbktmp = Cadd(Qbktmp, RCmul((double)((n + n + 1)*(1 - 2*(n % 2))), Csub(an[i], bn[i])));
-*/
-    //****************************************************//
-    // Calculate the scattering amplitudes (S1 and S2)    //
-    // Equations (25a) - (25b)                            //
-    //****************************************************//
-/*    for (t = 0; t < nTheta; t++) {
-      S1[t] = Cadd(S1[t], calc_S1(n, an[i], bn[i], Pi[i][t], Tau[i][t]));
-      S2[t] = Cadd(S2[t], calc_S2(n, an[i], bn[i], Pi[i][t], Tau[i][t]));
-    }*/
->>>>>>> eea51ce5ca0c5bb408c5a1c888b039637651b2e7
   }
 
-  return n_max;
+  return nmax;
 }
-

nmie.h

@@ -51,8 +51,8 @@ int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double>
          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 n_max,
-           int nCoords, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
-           std::vector<std::complex<double> > &E, std::vector<std::complex<double> >  &H);
+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);