merge

nmie.cc

@@ -363,8 +363,8 @@ namespace nmie {
   // Modify scattering (theta) angles                                       //
   // ********************************************************************** //
   void MultiLayerMie::SetAngles(const std::vector<double>& angles) {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     theta_ = angles;
   }
@@ -374,8 +374,8 @@ namespace nmie {
   // Modify size of all layers                                             //
   // ********************************************************************** //
   void MultiLayerMie::SetLayersSize(const std::vector<double>& layer_size) {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     size_param_.clear();
     double prev_layer_size = 0.0;
@@ -395,8 +395,8 @@ namespace nmie {
   // Modify refractive index of all layers                                  //
   // ********************************************************************** //
   void MultiLayerMie::SetLayersIndex(const std::vector< std::complex<double> >& index) {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     refractive_index_ = index;
   }
@@ -418,8 +418,8 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
   void MultiLayerMie::SetPECLayer(int layer_position) {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     if (layer_position < 0 && layer_position != -1)
       throw std::invalid_argument("Error! Layers are numbered from 0!");
@@ -431,8 +431,8 @@ namespace nmie {
   // Set maximun number of terms to be used                                 //
   // ********************************************************************** //
   void MultiLayerMie::SetMaxTerms(int nmax) {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     nmax_preset_ = nmax;
   }
@@ -453,8 +453,8 @@ namespace nmie {
   // Clear layer information                                                //
   // ********************************************************************** //
   void MultiLayerMie::ClearLayers() {
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
     size_param_.clear();
     refractive_index_.clear();
@@ -634,7 +634,7 @@ namespace nmie {
     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);
+      D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
     }
 
     if (std::abs(D1[0]) > 100000.0)
@@ -646,7 +646,7 @@ namespace nmie {
     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]);
+                                   *(static_cast<double>(n)*zinv - D3[n - 1]);
       D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
     }
     std::vector<std::complex<double> >  ZetaZ(nmax_+1), dZetaZ(nmax_ + 1);
@@ -684,7 +684,7 @@ namespace nmie {
     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]);
+      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]);
     }
     
@@ -712,9 +712,9 @@ namespace nmie {
   // using the relations:                                                             //
   //                                                                                  //
   //     j[n]   = Psi[n]/z                                                            //
-  //     j'[n]  = 0.5*(Psi[n-1]-Psi[n+1]-jn[n])/z                                     //
+  //     j'[n]  = j[n-1]-(n+1)*jn[n])/z                                               //
   //     h1[n]  = Zeta[n]/z                                                           //
-  //     h1'[n] = 0.5*(Zeta[n-1]-Zeta[n+1]-h1n[n])/z                                  //
+  //     h1'[n] = h1[n-1]-(n+1)*h1n[n]/z                                              //
   //                                                                                  //
   //**********************************************************************************//
   void MultiLayerMie::sbesjh(std::complex<double> z,
@@ -727,16 +727,16 @@ namespace nmie {
     calcPsiZeta(z, Psi, Zeta);
 
     // Now, calculate Spherical Bessel and Hankel functions and their derivatives
-    for (int n = 0; n < nmax_; n++) {
+    for (int n = 0; n <= nmax_; n++) {
       jn[n] = Psi[n]/z;
       h1n[n] = Zeta[n]/z;
 
       if (n == 0) {
-        jnp[n] = -Psi[1]/z - 0.5*jn[n]/z;
-        h1np[n] = -Zeta[1]/z - 0.5*h1n[n]/z;
+        jnp[0] = -Psi[1]/z - jn[0]/z;
+        h1np[0] = -Zeta[1]/z - h1n[0]/z;
       } else {
-        jnp[n] = 0.5*(Psi[n - 1] - Psi[n + 1] - jn[n])/z;
-        h1np[n] = 0.5*(Zeta[n - 1] - Zeta[n + 1] - h1n[n])/z;
+        jnp[n] = jn[n - 1] - static_cast<double>(n + 1)*jn[n]/z;
+        h1np[n] = h1n[n - 1] - static_cast<double>(n + 1)*h1n[n]/z;
       }
     }
   }
@@ -793,7 +793,7 @@ namespace nmie {
   // Output parameters:                                                               //
   //   Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics                     //
   //**********************************************************************************//
-  void MultiLayerMie::calcSpherHarm(const double Rho, const double Phi, const double Theta,
+  void MultiLayerMie::calcSpherHarm(const double Rho, const double Theta, const double Phi,
                                     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, 
@@ -839,9 +839,9 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  void MultiLayerMie::ExternalScattCoeffs() {
+  void MultiLayerMie::ScattCoeffs() {
 
-    isExtCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
 
     const std::vector<double>& x = size_param_;
     const std::vector<std::complex<double> >& m = refractive_index_;
@@ -988,8 +988,8 @@ namespace nmie {
         bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
       }
     }  // end of for an and bn terms
-    isExtCoeffsCalc_ = true;
-  }  // end of MultiLayerMie::ExternalScattCoeffs(...)
+    isScaCoeffsCalc_ = true;
+  }  // end of MultiLayerMie::ScattCoeffs(...)
 
 
   //**********************************************************************************//
@@ -1028,14 +1028,14 @@ namespace nmie {
 
     const std::vector<double>& x = size_param_;
 
-    isIntCoeffsCalc_ = false;
-    isExtCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
+    isScaCoeffsCalc_ = false;
     isMieCalculated_ = false;
 
     // Calculate scattering coefficients
-    ExternalScattCoeffs();
+    ScattCoeffs();
 
-    if (!isExtCoeffsCalc_) // TODO seems to be unreachable
+    if (!isScaCoeffsCalc_) // TODO seems to be unreachable
       throw std::invalid_argument("Calculation of scattering coefficients failed!");
 
     // Initialize the scattering parameters
@@ -1094,7 +1094,7 @@ namespace nmie {
 
 
   //**********************************************************************************//
-  // This function calculates the scattering coefficients inside the particle,        //
+  // This function calculates the expansion coefficients inside the particle,         //
   // required to calculate the near-field parameters.                                 //
   //                                                                                  //
   // Input parameters:                                                                //
@@ -1112,14 +1112,13 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  void MultiLayerMie::InternalScattCoeffs() {
-    if (!isExtCoeffsCalc_)
-      throw std::invalid_argument("(InternalScattCoeffs) You should calculate external coefficients first!");
+  void MultiLayerMie::ExpanCoeffs() {
+    if (!isScaCoeffsCalc_)
+      throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
 
-    isIntCoeffsCalc_ = false;
+    isExpCoeffsCalc_ = false;
 
-    std::complex<double> c_one(1.0, 0.0);
-    std::complex<double> c_zero(0.0, 0.0);
+    std::complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
 
     const int L = refractive_index_.size();
 
@@ -1127,72 +1126,62 @@ namespace nmie {
     bln_.resize(L + 1);
     cln_.resize(L + 1);
     dln_.resize(L + 1);
-    for (auto& element:aln_) element.resize(nmax_);
-    for (auto& element:bln_) element.resize(nmax_);
-    for (auto& element:cln_) element.resize(nmax_);
-    for (auto& element:dln_) element.resize(nmax_);
+    for (int l = 0; l <= L; l++) {
+      aln_[l].resize(nmax_);
+      bln_[l].resize(nmax_);
+      cln_[l].resize(nmax_);
+      dln_[l].resize(nmax_);
+    }
 
     // Yang, paragraph under eq. A3
     // a^(L + 1)_n = a_n, d^(L + 1) = 1 ...
-    for (int i = 0; i < nmax_; ++i) {
-      aln_[L][i] = an_[i];
-      bln_[L][i] = bn_[i];
-      cln_[L][i] = c_one;
-      dln_[L][i] = c_one;
+    for (int n = 0; n < nmax_; n++) {
+      aln_[L][n] = an_[n];
+      bln_[L][n] = bn_[n];
+      cln_[L][n] = c_one;
+      dln_[L][n] = c_one;
     }
 
-    std::vector<std::complex<double> > z(L), z1(L);
-    for (int i = 0; i < L - 1; ++i) {
-      z[i] = size_param_[i]*refractive_index_[i];
-      z1[i] = size_param_[i]*refractive_index_[i + 1];
-    }
-    z[L - 1] = size_param_[L - 1]*refractive_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);
-    }
+    std::vector<std::complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
+    std::vector<std::complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
+    std::complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
 
-    for (int l = 0; l < L; ++l) {
-      calcD1D3(z[l], D1z[l], D3z[l]);
-      calcD1D3(z1[l], D1z1[l], D3z1[l]);
-      calcPsiZeta(z[l], Psiz[l], Zetaz[l]);
-      calcPsiZeta(z1[l], Psiz1[l], Zetaz1[l]);
-    }
     auto& m = refractive_index_;
     std::vector< std::complex<double> > m1(L);
 
-    for (int l = 0; l < L - 1; ++l) m1[l] = m[l + 1];
+    for (int l = 0; l < L - 1; l++) m1[l] = m[l + 1];
     m1[L - 1] = std::complex<double> (1.0, 0.0);
 
+    std::complex<double> z, z1;
     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 = aln_[l + 1][n]*Zetaz1[l][n + 1] - dln_[l + 1][n]*Psiz1[l][n + 1];
-        auto T2 = bln_[l + 1][n]*Zetaz1[l][n + 1] - cln_[l + 1][n]*Psiz1[l][n + 1];
-
-        auto T3 = -D1z1[l][n + 1]*dln_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*aln_[l + 1][n]*Zetaz1[l][n + 1];
-        auto T4 = -D1z1[l][n + 1]*cln_[l + 1][n]*Psiz1[l][n + 1] + D3z1[l][n + 1]*bln_[l + 1][n]*Zetaz1[l][n + 1];
-
-        // anl
-        aln_[l][n] = (D1z[l][n + 1]*m1[l]*T1 - m[l]*T3)/denomZeta;
-        // bnl
-        bln_[l][n] = (D1z[l][n + 1]*m[l]*T2 - m1[l]*T4)/denomZeta;
-        // cnl
-        cln_[l][n] = (D3z[l][n + 1]*m[l]*T2 - m1[l]*T4)/denomPsi;
-        // dnl
-        dln_[l][n] = (D3z[l][n + 1]*m1[l]*T1 - m[l]*T3)/denomPsi;
+      z = size_param_[l]*m[l];
+      z1 = size_param_[l]*m1[l];
+
+      calcD1D3(z, D1z, D3z);
+      calcD1D3(z1, D1z1, D3z1);
+      calcPsiZeta(z, Psiz, Zetaz);
+      calcPsiZeta(z1, Psiz1, Zetaz1);
+
+      for (int n = 0; n < nmax_; n++) {
+        int n1 = n + 1;
+
+        denomZeta = m1[l]*Zetaz[n1]*(D1z[n1] - D3z[n1]);
+        denomPsi  =  m1[l]*Psiz[n1]*(D1z[n1] - D3z[n1]);
+
+        T1 = aln_[l + 1][n]*Zetaz1[n1] - dln_[l + 1][n]*Psiz1[n1];
+        T2 = bln_[l + 1][n]*Zetaz1[n1] - cln_[l + 1][n]*Psiz1[n1];
+
+        T3 = D1z1[n1]*dln_[l + 1][n]*Psiz1[n1] - D3z1[n1]*aln_[l + 1][n]*Zetaz1[n1];
+        T4 = D1z1[n1]*cln_[l + 1][n]*Psiz1[n1] - D3z1[n1]*bln_[l + 1][n]*Zetaz1[n1];
+
+        // aln
+        aln_[l][n] = (D1z[n1]*m1[l]*T1 + m[l]*T3)/denomZeta;
+        // bln
+        bln_[l][n] = (D1z[n1]*m[l]*T2 + m1[l]*T4)/denomZeta;
+        // cln
+        cln_[l][n] = (D3z[n1]*m[l]*T2 + m1[l]*T4)/denomPsi;
+        // dln
+        dln_[l][n] = (D3z[n1]*m1[l]*T1 + m[l]*T3)/denomPsi;
       }  // end of all n
     }  // end of all l
 
@@ -1206,8 +1195,8 @@ namespace nmie {
       else throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
     }
 
-    isIntCoeffsCalc_ = true;
-  }  // end of   void MultiLayerMie::InternalScattCoeffs()
+    isExpCoeffsCalc_ = true;
+  }  // end of   void MultiLayerMie::ExpanCoeffs()
 
 
   // ********************************************************************** //
@@ -1215,7 +1204,7 @@ namespace nmie {
   // BH p.92 (4.37), 94 (4.45), 95 (4.50)                                   //
   // assume: medium is non-absorbing; refim = 0; Uabs = 0                   //
   // ********************************************************************** //
-  void MultiLayerMie::fieldExt(const double Rho, const double Phi, const double Theta,
+  void MultiLayerMie::fieldExt(const double Rho, const double Theta, const double Phi,
                                std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 
     std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
@@ -1236,7 +1225,7 @@ namespace nmie {
       double rn = static_cast<double>(n1);
 
       // using BH 4.12 and 4.50
-      calcSpherHarm(Rho, Phi, Theta, h1n[n1], h1np[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
+      calcSpherHarm(Rho, Theta, Phi, h1n[n1], h1np[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
 
       // scattered field: BH p.94 (4.45)
       std::complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
@@ -1294,20 +1283,25 @@ namespace nmie {
   // Output parameters:                                                               //
   //   E, H: Complex electric and magnetic fields                                     //
   //**********************************************************************************//
-  void MultiLayerMie::calcField(const double Rho, const double Phi, const double Theta,
+  void MultiLayerMie::calcField(const double Rho, const double Theta, const double Phi,
                                 std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 
     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), 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 l = 0;  // Layer number
     std::complex<double> ml;
 
+    // Initialize E and H
+    for (int i = 0; i < 3; i++) {
+      E[i] = c_zero;
+      H[i] = c_zero;
+    }
+
     if (Rho > size_param_.back()) {
       l = size_param_.size();
       ml = c_one;
@@ -1326,35 +1320,32 @@ namespace nmie {
     // Calculate angular functions Pi and Tau
     calcPiTau(std::cos(Theta), Pi, Tau);
 
-    for (int n = nmax_ - 1; n >= 0; n--) {
+    for (int n = nmax_ - 2; n >= 0; n--) {
       int n1 = n + 1;
       double rn = static_cast<double>(n1);
 
       // 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[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
+      calcSpherHarm(Rho, Theta, Phi, jn[n1], jnp[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
+      calcSpherHarm(Rho, Theta, Phi, h1n[n1], h1np[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
 
       // 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++) {
-        El[i] = El[i] + En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
-                      + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
+        // electric field E [V m - 1] = EF*E0
+        E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
+              + c_i*aln_[l][n]*N3e1n[i] -     bln_[l][n]*M3o1n[i]);
 
-        Hl[i] = Hl[i] + En*(-dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
-                      +  c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
+        H[i] += En*(-dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
+              +  c_i*bln_[l][n]*N3o1n[i] +     aln_[l][n]*M3e1n[i]);
       }
     }  // end of for all n
 
+    //printf("rho = %10.5e; phi = %10.5eº; theta = %10.5eº; x[%i] = %10.5e; m[%i] = %10.5er%+10.5ei\n", Rho, Phi*180./PI_, Theta*180./PI_, l, size_param_[l], l, std::real(ml), std::imag(ml));
     // magnetic field
     double hffact = 1.0/(cc_*mu_);
     for (int i = 0; i < 3; i++) {
-      Hl[i] = hffact*Hl[i];
-    }
-
-    for (int i = 0; i < 3; i++) {
-      // electric field E [V m - 1] = EF*E0
-      E[i] = El[i];
-      H[i] = Hl[i];
+      H[i] = hffact*H[i];
+      //printf("E[%i] = %10.5er%+10.5ei; H[%i] = %10.5er%+10.5ei\n", i, std::real(E[i]), std::imag(E[i]), i, std::real(H[i]), std::imag(H[i]));
     }
    }  // end of MultiLayerMie::calcField(...)
 
@@ -1382,22 +1373,20 @@ namespace nmie {
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
   void MultiLayerMie::RunFieldCalculation() {
-    // Calculate external scattering coefficients an_ and bn_
-    ExternalScattCoeffs();
-    printf("an_\n");
-    for (auto a: an_) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
-    printf("\nbn_\n");
-    for (auto a: bn_) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
-
-    calc_an_bn_bulk(an_, bn_, size_param_.back(), refractive_index_.back());
-    printf("\nbulk an_\n");
-    for (auto a: an_) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
-    printf("\nbulk bn_\n");
-    for (auto a: bn_) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
-    printf("\nsize param: %g\n", size_param_.back());
-
-    // Calculate internal scattering coefficients aln_,  bln_, cln_, and dln_
-    InternalScattCoeffs();
+    double Rho, Theta, Phi;
+
+    // Calculate scattering coefficients an_ and bn_
+    ScattCoeffs();
+
+    std::vector<std::complex<double> > an1(nmax_), bn1(nmax_);
+    calc_an_bn_bulk(an1, bn1, size_param_.back(), refractive_index_.back());
+    for (int n = 0; n < nmax_; n++) {
+      printf("an_[%i] = %10.5er%+10.5ei;  an_bulk_[%i] = %10.5er%+10.5ei\n", n, std::real(an_[n]), std::imag(an_[n]), n, std::real(an1[n]), std::imag(an1[n]));
+      printf("bn_[%i] = %10.5er%+10.5ei;  bn_bulk_[%i] = %10.5er%+10.5ei\n", n, std::real(bn_[n]), std::imag(bn_[n]), n, std::real(bn1[n]), std::imag(bn1[n]));
+    }
+
+    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
+    ExpanCoeffs();
 
     long total_points = coords_[0].size();
     E_.resize(total_points);
@@ -1411,17 +1400,22 @@ namespace nmie {
       const double& Zp = coords_[2][point];
 
       // Convert to spherical coordinates
-      double Rho = std::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
+      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;
+      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;
+      if (Xp == 0.0)
+        Phi = (Yp != 0.0) ? std::asin(Yp/std::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
+      else
+        Phi = std::acos(Xp/std::sqrt(pow2(Xp) + pow2(Yp)));
 
       // Avoid convergence problems due to Rho too small
       if (Rho < 1e-5) Rho = 1e-5;
 
+      //printf("X = %g; Y = %g; Z = %g; pho = %g; phi = %g; theta = %g\n", Xp, Yp, Zp, Rho, Phi*180./PI_, Theta*180./PI_);
+
       //*******************************************************//
       // external scattering field = incident + scattered      //
       // BH p.92 (4.37), 94 (4.45), 95 (4.50)                  //
@@ -1433,9 +1427,9 @@ namespace nmie {
 
       // Firstly the easiest case: the field outside the particle
       //      if (Rho >= GetSizeParameter()) {
-      //        fieldExt(Rho, Phi, Theta, Es, Hs);
+      //        fieldExt(Rho, Theta, Phi, Es, Hs);
       // } else {
-      calcField(Rho, Phi, Theta, Es, Hs);  //Should work fine both: inside and outside the particle
+      calcField(Rho, Theta, Phi, Es, Hs);  //Should work fine both: inside and outside the particle
       //}
 
       { //Now, convert the fields back to cartesian coordinates

nmie.h

@@ -124,22 +124,22 @@ namespace nmie {
                 std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np);
     void calcPiTau(const double& costheta,
                    std::vector<double>& Pi, std::vector<double>& Tau);
-    void calcSpherHarm(const double Rho, const double Phi, const double Theta,
+    void calcSpherHarm(const double Rho, const double Theta, const double Phi,
                        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 ExternalScattCoeffs();
-    void InternalScattCoeffs();
+    void ScattCoeffs();
+    void ExpanCoeffs();
 
-    void fieldExt(const double Rho, const double Phi, const double Theta,
+    void fieldExt(const double Rho, const double Theta, const double Phi,
                   std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
 
-    void calcField(const double Rho, const double Phi, const double Theta,
+    void calcField(const double Rho, const double Theta, const double Phi,
                    std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
 
-    bool isIntCoeffsCalc_ = false;
-    bool isExtCoeffsCalc_ = false;
+    bool isExpCoeffsCalc_ = false;
+    bool isScaCoeffsCalc_ = false;
     bool isMieCalculated_ = false;
 
     // Size parameter for all layers

scattnlay.cpp

@@ -1,4 +1,4 @@
-/* Generated by Cython 0.20.1post0 (Debian 0.20.1+git90-g0e6e38e-1ubuntu2) on Tue Apr 14 15:24:18 2015 */
+/* Generated by Cython 0.20.1post0 (Debian 0.20.1+git90-g0e6e38e-1ubuntu2) on Wed Apr 15 19:02:11 2015 */
 
 #define PY_SSIZE_T_CLEAN
 #ifndef CYTHON_USE_PYLONG_INTERNALS

setup.py

@@ -53,7 +53,7 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
       license = 'GPL',
       platforms = 'any',
       ext_modules = [Extension("scattnlay",
-                               ["nmie.cc", "bessel.cc","py_nmie.cc", "scattnlay.cpp"],
+                               ["nmie.cc", "bessel.cc", "py_nmie.cc", "scattnlay.cpp"],
                                language = "c++",
                                include_dirs = [np.get_include()])], 
       extra_compile_args=['-std=c++11']

tests/python/field.py

@@ -35,9 +35,9 @@
 # conditions.
 import scattnlay
 
-import os
-path = os.path.dirname(scattnlay.__file__)
-print(scattnlay.__file__)
+#import os
+#path = os.path.dirname(scattnlay.__file__)
+#print(scattnlay.__file__)
 
 from scattnlay import fieldnlay
 import numpy as np
@@ -52,7 +52,7 @@ m[0, 1] = (0.565838 + 7.23262j)/1.3205
 
 npts = 501
 
-scan = np.linspace(-2.0*x[0, 1], 2.0*x[0, 1], 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)

tests/python/lfield.py

@@ -58,12 +58,6 @@ 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

@@ -56,7 +56,7 @@ 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 "x =", x
+print "m =", m
 print np.vstack((coord[:, 0], Eh)).transpose()