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

nmie-wrapper.cc

@@ -72,7 +72,7 @@ namespace nmie {
       S1 = multi_layer_mie.GetS1();
       S2 = multi_layer_mie.GetS2();
       //multi_layer_mie.GetFailed();
-    } catch( const std::invalid_argument& ia ) {
+    } 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);
@@ -86,16 +86,16 @@ namespace nmie {
   // ********************************************************************** //
   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!");
+      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 )
+        || 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!");
+      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!");
+      if (f.size() != 3)
+        throw std::invalid_argument("Field H is not 3D!");
     try {
       MultiLayerMie multi_layer_mie;  
       //multi_layer_mie.SetPEC(pl);
@@ -106,11 +106,11 @@ namespace nmie {
       E = multi_layer_mie.GetFieldE();
       H = multi_layer_mie.GetFieldH();
       //multi_layer_mie.GetFailed();
-    } catch( const std::invalid_argument& ia ) {
+    } 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 - 1;
     }  
 
     return 0;
@@ -124,18 +124,18 @@ namespace nmie {
     std::complex<double> z(faild_x, 0.0);
     std::vector<int> nmax_local_array = {20, 100, 500, 2500};
     for (auto nmax_local : nmax_local_array) {
-      std::vector<std::complex<double> > D1_failed(nmax_local +1);
+      std::vector<std::complex<double> > D1_failed(nmax_local + 1);
       // Downward recurrence for D1 - equations (16a) and (16b)
       D1_failed[nmax_local] = std::complex<double>(0.0, 0.0);
       const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
       for (int n = nmax_local; n > 0; n--) {
-	D1_failed[n - 1] = double(n)*zinv - 1.0/(D1_failed[n] + double(n)*zinv);
+        D1_failed[n - 1] = double(n)*zinv - 1.0/(D1_failed[n] + double(n)*zinv);
       }
       printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n",
-	     faild_x, nmax_local, D1_failed[0].real());
+             faild_x, nmax_local, D1_failed[0].real());
     }
     printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x,
-	   calcD1confra(0,z).real());
+           calcD1confra(0,z).real());
     //D1[nmax_] = calcD1confra(nmax_, z);
   
     
@@ -170,13 +170,13 @@ namespace nmie {
   std::vector<double> MultiLayerMie::GetQabs_channel_normalized() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    // std::vector<double> NACS(nmax_-1, 0.0);
+    // std::vector<double> NACS(nmax_ - 1, 0.0);
     // double x2 = pow2(size_parameter_.back());
     // for (int i = 0; i < nmax_ - 1; ++i) {
-    //   const int n = i+1;
-    //   NACS[i] = Qabs_ch_[i]*x2/(2.0*(2.0*static_cast<double>(n)+1));
+    //   const int n = i + 1;
+    //   NACS[i] = Qabs_ch_[i]*x2/(2.0*(2.0*static_cast<double>(n) + 1));
     //   // if (NACS[i] > 0.250000001)
-    //   // 	throw std::invalid_argument("Unexpected normalized absorption cross-section value!");
+    //   //         throw std::invalid_argument("Unexpected normalized absorption cross-section value!");
     // }
     //return NACS;    
     return Qabs_ch_norm_;
@@ -203,11 +203,11 @@ namespace nmie {
   std::vector<double> MultiLayerMie::GetQsca_channel_normalized() {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
-    // std::vector<double> NACS(nmax_-1, 0.0);
+    // std::vector<double> NACS(nmax_ - 1, 0.0);
     // double x2 = pow2(size_parameter_.back());
     // for (int i = 0; i < nmax_ - 1; ++i) {
-    //   const int n = i+1;
-    //   NACS[i] = Qsca_ch_[i]*x2/(2.0*(2.0*static_cast<double>(n)+1.0));
+    //   const int n = i + 1;
+    //   NACS[i] = Qsca_ch_[i]*x2/(2.0*(2.0*static_cast<double>(n) + 1.0));
     // }
     // return NACS;    
      return Qsca_ch_norm_;
@@ -323,7 +323,7 @@ namespace nmie {
         throw std::invalid_argument("Size parameter should be positive!");
       if (prev_size_parameter > layer_size_parameter) 
         throw std::invalid_argument
-	  ("Size parameter for next layer should be larger than the previous one!");
+          ("Size parameter for next layer should be larger than the previous one!");
       prev_size_parameter = layer_size_parameter;
       size_parameter_.push_back(layer_size_parameter);
     }
@@ -344,8 +344,7 @@ namespace nmie {
   void MultiLayerMie::SetFieldPointsSP(const std::vector< std::vector<double> >& coords_sp) {
     if (coords_sp.size() != 3)
       throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
-    if (coords_sp[0].size() != coords_sp[1].size()
-	|| coords_sp[0].size() != coords_sp[2].size())
+    if (coords_sp[0].size() != coords_sp[1].size() || coords_sp[0].size() != coords_sp[2].size())
       throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
     coords_sp_ = coords_sp;
     // for (int i = 0; i < coords_sp_[0].size(); ++i) {
@@ -379,11 +378,11 @@ namespace nmie {
     double radius = 0.0;
     for (auto width : target_width_) {
       radius += width;
-      size_parameter_.push_back(2*PI_*radius / wavelength_);
+      size_parameter_.push_back(2*PI_*radius/wavelength_);
     }
     for (auto width : coating_width_) {
       radius += width;
-      size_parameter_.push_back(2*PI_*radius / wavelength_);
+      size_parameter_.push_back(2*PI_*radius/wavelength_);
     }
     total_radius_ = radius;
   }  // end of void MultiLayerMie::GenerateSizeParameter();
@@ -413,14 +412,14 @@ namespace nmie {
     if (!isMieCalculated_)
       throw std::invalid_argument("You should run calculations before result request!");
     std::vector< std::vector<double> > spectra;
-    double step_WL = (to_WL - from_WL)/ static_cast<double>(samples);
+    double step_WL = (to_WL - from_WL)/static_cast<double>(samples);
     double wavelength_backup = wavelength_;
     long fails = 0;
     for (double WL = from_WL; WL < to_WL; WL += step_WL) {
       wavelength_ = WL;
       try {
         RunMieCalculations();
-      } catch( const std::invalid_argument& ia ) {
+      } catch(const std::invalid_argument& ia) {
         fails++;
         continue;
       }
@@ -493,15 +492,15 @@ namespace nmie {
     Nstop();  // 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]));
+        ri = round(std::abs(x[i]*m[i]));
       else 
-	ri = 0;      
+        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]));
+        riM1 = round(std::abs(x[i - 1]* m[i]));
       else 
-	riM1 = 0;      
+        riM1 = 0;      
       nmax_ = std::max(nmax_, riM1);
     }
     nmax_ += 15;  // Final nmax_ value
@@ -521,7 +520,7 @@ namespace nmie {
   // 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      //
+  // 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               //
@@ -533,10 +532,10 @@ namespace nmie {
   // 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) {
+                             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;
@@ -565,12 +564,12 @@ namespace nmie {
 
       absc = std::abs(std::real(d)) + std::abs(std::imag(d));
       if (absc < tm30) {
-	d = tm30;
+        d = tm30;
       }
 
       absc = std::abs(std::real(c)) + std::abs(std::imag(c));
       if (absc < tm30) {
-	c = tm30;
+        c = tm30;
       }
 
       d = 1.0/d;
@@ -581,8 +580,8 @@ namespace nmie {
       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;
+        // We have obtained the desired accuracy
+        break;
       }
     }
 
@@ -613,18 +612,18 @@ namespace nmie {
       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];
+        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];
-	}
+        // 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;
+        pl += zi;
 
-	h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
+        h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
       }
     }
   }
@@ -642,9 +641,9 @@ namespace nmie {
   //   bd: Logarithmic derivative                                                     //
   //**********************************************************************************//
   void MultiLayerMie::sphericalBessel(std::complex<double> z,
-				      std::vector<std::complex<double> >& bj,
-				      std::vector<std::complex<double> >& by,
-				      std::vector<std::complex<double> >& bd) {
+                                      std::vector<std::complex<double> >& bj,
+                                      std::vector<std::complex<double> >& by,
+                                      std::vector<std::complex<double> >& bd) {
     std::vector<std::complex<double> > jn(nmax_), jnp(nmax_), h1n(nmax_), h1np(nmax_);
     sbesjh(z, jn, jnp, h1n, h1np);
 
@@ -664,8 +663,8 @@ namespace nmie {
     //   by[1] = by[0]*3.0/z-bessely_0;//bj2
     // }
     // for (int n = 2; n < nmax_; n++) {
-    //   bj[n] = (2.0*n-1.0)/z*bj[n-1] - bj[n];
-    //   by[n] = (2.0*n-1.0)/z*by[n-1] - by[n];
+    //   bj[n] = (2.0*n - 1.0)/z*bj[n - 1] - bj[n];
+    //   by[n] = (2.0*n - 1.0)/z*by[n - 1] - by[n];
     // }
   }
   // ********************************************************************** //
@@ -673,8 +672,8 @@ namespace nmie {
   // ********************************************************************** //
   // 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::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;
@@ -686,8 +685,8 @@ namespace nmie {
   // ********************************************************************** //
   // 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> 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;
@@ -699,7 +698,7 @@ namespace nmie {
   // ********************************************************************** //
   // Calculates S1 - equation (25a)
   std::complex<double> MultiLayerMie::calc_S1(int n, std::complex<double> an, std::complex<double> bn,
-					      double Pi, double Tau) {
+                                              double Pi, double Tau) {
     return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
   }
   // ********************************************************************** //
@@ -707,7 +706,7 @@ namespace nmie {
   // ********************************************************************** //
   // 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) {
+                                              double Pi, double Tau) {
     return calc_S1(n, an, bn, Tau, Pi);
   }
   //**********************************************************************************//
@@ -723,10 +722,10 @@ namespace nmie {
   //   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) {
+                                  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)
     //Psi[0] = std::complex<double>(std::sin(x), 0);
     std::complex<double> c_i(0.0, 1.0);
@@ -758,33 +757,33 @@ c     a_k
 c    CAPT     Factor used in Lentz iteration for A (Eq. R27)
 c     T_k
 
-c    CNUMER   Numerator   in capT  ( Eq. R28A )
+c    CNUMER   Numerator   in capT  (Eq. R28A)
 c     N_k
-c    CDENOM   Denominator in capT  ( Eq. R28B )
+c    CDENOM   Denominator in capT  (Eq. R28B)
 c     D_k
 
 c    CDTD     Product of two successive denominators of capT factors
-c                 ( Eq. R34C )
+c                 (Eq. R34C)
 c     xi_1
 
 c    CNTN     Product of two successive numerators of capT factors
-c                 ( Eq. R34B )
+c                 (Eq. R34B)
 c     xi_2
 
 c    EPS1     Ill-conditioning criterion
 c    EPS2     Convergence criterion
 
-c    KK       Subscript k of cAk  ( Eq. R25B )
+c    KK       Subscript k of cAk  (Eq. R25B)
 c     k
 
-c    KOUNT    Iteration counter ( used to prevent infinite looping )
+c    KOUNT    Iteration counter (used to prevent infinite looping)
 
 c    MAXIT    Max. allowed no. of iterations
 
-c    MM       + 1  and - 1, alternately
+c    MM + 1  and - 1, alternately
 */
   std::complex<double> MultiLayerMie::calcD1confra(const int N, const std::complex<double> z) {
-  // NTMR -> nmax_ -1  \\TODO nmax_ ?
+  // NTMR -> nmax_ - 1  \\TODO nmax_ ?
     //int N = nmax_ - 1;
     int KK, KOUNT, MAXIT = 10000, MM;
     //    double EPS1=1.0e-2;
@@ -793,49 +792,49 @@ c    MM       + 1  and - 1, alternately
     std::complex<double> one = std::complex<double>(1.0,0.0);
     std::complex<double> ZINV = one/z;
 // c                                 ** Eq. R25a
-    std::complex<double> CONFRA = static_cast<std::complex<double> >(N+1)*ZINV;   //debug ZINV
-    MM = -1; 
+    std::complex<double> CONFRA = static_cast<std::complex<double> >(N + 1)*ZINV;   //debug ZINV
+    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
     CDENOM = CAK;
-    CNUMER = CDENOM + one / CONFRA; //-3zinv+z
+    CNUMER = CDENOM + one/CONFRA; //-3zinv+z
     KOUNT  = 1;
     //10 CONTINUE
     do {      ++KOUNT;
       if (KOUNT > MAXIT) {
-	printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
-	throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
+        printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
+        throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
       }
-      MM *= -1;      KK += 2;  //debug  mm=1 kk=5
+      MM *= - 1;      KK += 2;  //debug  mm=1 kk=5
       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
-     // 	MM  *= -1;	KK  += 2;
-     // 	CAK = static_cast<std::complex<double> >(MM*KK) * ZINV; // ** Eq. R25b
-     // 	//c                        ** Eq. R35
-     // 	CNUMER = CAK + CNUMER / CNTN;
-     // 	CDENOM = CAK + CDENOM / CDTD;
-     // 	++KOUNT;
-     // 	//GO TO  10
-     // 	continue;
+     //  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
+     //         MM  *= - 1;        KK  += 2;
+     //         CAK = static_cast<std::complex<double> >(MM*KK) * ZINV; // ** Eq. R25b
+     //         //c                        ** Eq. R35
+     //         CNUMER = CAK + CNUMER/CNTN;
+     //         CDENOM = CAK + CDENOM/CDTD;
+     //         ++KOUNT;
+     //         //GO TO  10
+     //         continue;
      // } else { //c                           *** Well-conditioned case
       {
-	CAPT   = CNUMER / CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
-	// printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
+        CAPT   = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
+        // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
        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 ) {
+       if (std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2) {
 //c                                        ** Eq. R30
-	 CNUMER = CAK + one/CNUMER;
-	 CDENOM = CAK + one/CDENOM;
-	 continue;
-	 //GO TO  10
+         CNUMER = CAK + one/CNUMER;
+         CDENOM = CAK + one/CDENOM;
+         continue;
+         //GO TO  10
        }  // end of if < eps2
       }
       break;
@@ -856,8 +855,8 @@ c    MM       + 1  and - 1, alternately
   //   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) {
+                               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);
     //D1[nmax_] = calcD1confra(nmax_, z);
@@ -867,21 +866,21 @@ c    MM       + 1  and - 1, alternately
     // prn((D1[nmax_] + double(nmax_)*zinv).real());
     for (int n = nmax_; n > 0; n--) {
       D1[n - 1] = double(n)*zinv - 1.0/(D1[n] + double(n)*zinv);
-      //D1[n-1] = calcD1confra(n-1, z);
-      // printf(" D:");prn((D1[n-1]).real()); printf("\t diff:");
+      //D1[n - 1] = calcD1confra(n - 1, z);
+      // printf(" D:");prn((D1[n - 1]).real()); printf("\t diff:");
       // prn((D1[n] + double(n)*zinv).real());
     }
     //     printf("\n\n"); iformat=0;
-    if (std::abs(D1[0]) > 100000.0 )
+    if (std::abs(D1[0]) > 100000.0)
       throw std::invalid_argument
-	("Unstable D1! Please, try to change input parameters!\n");
+        ("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()));
+                       *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]);
+        *(static_cast<double>(n)*zinv- D3[n - 1]);
       D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
     }
   }
@@ -890,7 +889,7 @@ c    MM       + 1  and - 1, alternately
   // Equations (26a) - (26c)                                                          //
   //                                                                                  //
   // Input parameters:                                                                //
-  //   nmax_: 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                                           //
@@ -899,35 +898,35 @@ c    MM       + 1  and - 1, alternately
   //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
   //**********************************************************************************//
   void MultiLayerMie::calcSinglePiTau(const double& costheta, std::vector<double>& Pi,
-				      std::vector<double>& Tau) {
+                                      std::vector<double>& Tau) {
     //****************************************************//
     // Equations (26a) - (26c)                            //
     //****************************************************//
     for (int n = 0; n < nmax_; n++) {
       if (n == 0) {
-	// Initialize Pi and Tau
-	Pi[n] = 1.0;
-	Tau[n] = (n + 1)*costheta; 
+        // Initialize Pi and Tau
+        Pi[n] = 1.0;
+        Tau[n] = (n + 1)*costheta; 
       } else {
-	// Calculate the actual values
-	Pi[n] = ((n == 1) ? ((n + n + 1)*costheta*Pi[n - 1]/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];
+        // Calculate the actual values
+        Pi[n] = ((n == 1) ? ((n + n + 1)*costheta*Pi[n - 1]/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 void MultiLayerMie::calcPiTau(...)
   void MultiLayerMie::calcAllPiTau(std::vector< std::vector<double> >& Pi,
-				std::vector< std::vector<double> >& Tau) {
+                                   std::vector< std::vector<double> >& Tau) {
     std::vector<double> costheta(theta_.size(), 0.0);
-    for (int t = 0; t < theta_.size(); t++) {	
+    for (int t = 0; t < theta_.size(); t++) {
       costheta[t] = std::cos(theta_[t]);
     }
     // Do not join upper and lower 'for' to a single one!  It will slow
     // down the code!!! (For about 0.5-2.0% of runtime, it is probably
     // due to increased cache missing rate originated from the
     // recurrence in calcPiTau...)
-    for (int t = 0; t < theta_.size(); t++) {	
+    for (int t = 0; t < theta_.size(); t++) {
       calcSinglePiTau(costheta[t], Pi[t], Tau[t]);
       //calcSinglePiTau(std::cos(theta_[t]), Pi[t], Tau[t]); // It is slow!!
     }
@@ -952,19 +951,19 @@ c    MM       + 1  and - 1, alternately
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
   void MultiLayerMie::ScattCoeffs(std::vector<std::complex<double> >& an,
-				  std::vector<std::complex<double> >& bn) {
+                                  std::vector<std::complex<double> >& bn) {
     const std::vector<double>& x = size_parameter_;
     const std::vector<std::complex<double> >& m = index_;
     const int& pl = PEC_layer_position_;
     const int L = 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.  //
-    // ************************************************************************//
+    // 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;
+    // 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.
     int fl = (pl > 0) ? pl : 0;
@@ -1001,8 +1000,8 @@ c    MM       + 1  and - 1, alternately
     //*************************************************//
     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);
+        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];
@@ -1011,7 +1010,7 @@ c    MM       + 1  and - 1, alternately
     }
     // do { \
     //   ++iformat;\
-    //   if (iformat%5 == 0) printf("%24.16e",z1.real());	\
+    //   if (iformat%5 == 0) printf("%24.16e",z1.real());
     // } while (false);
     //******************************************************************//
     // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
@@ -1027,7 +1026,7 @@ c    MM       + 1  and - 1, alternately
     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     //
+      //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];
@@ -1036,50 +1035,50 @@ c    MM       + 1  and - 1, alternately
       //Calculate D1 and D3 for z2
       calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
       // prn(z1.real());
-      // for ( auto i : D1_mlxl) { prn(i.real());
+      // for (auto i : D1_mlxl) { prn(i.real());
       //   // prn(i.imag());
-      // 	} printf("\n");
+      //         } printf("\n");
 
       //*********************************************//
-      //Calculate Q, Ha and Hb in the layers fl+1..L //
+      //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()));
+        * 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;
+        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
+        //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);
+        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
     //**************************************//
@@ -1087,7 +1086,7 @@ c    MM       + 1  and - 1, alternately
     //**************************************//
     // Calculate D1XL and D3XL
     calcD1D3(x[L - 1],  D1XL, D3XL);
-    //printf("%5.20f\n",Ha[L-1][0].real());
+    //printf("%5.20f\n",Ha[L - 1][0].real());
     // Calculate PsiXL and ZetaXL
     calcPsiZeta(x[L - 1], D1XL, D3XL, PsiXL, ZetaXL);
     //*********************************************************************//
@@ -1102,11 +1101,11 @@ c    MM       + 1  and - 1, alternately
       //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]);
+        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];
+        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
   }  // end of void MultiLayerMie::ScattCoeffs(...)
@@ -1128,18 +1127,18 @@ c    MM       + 1  and - 1, alternately
     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);
+    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));
+    std::vector<std::complex<double> >        tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
     S1_.swap(tmp1);
     S2_ = S1_;
   }
@@ -1163,15 +1162,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          //
   //         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                                         //
@@ -1222,11 +1221,11 @@ c    MM       + 1  and - 1, alternately
       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());
+                          + 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());
+      //                             + 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
@@ -1234,7 +1233,7 @@ c    MM       + 1  and - 1, alternately
       // 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());
+               + ((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]);
@@ -1242,8 +1241,8 @@ c    MM       + 1  and - 1, alternately
       // Calculate the scattering amplitudes (S1 and S2)    //
       // Equations (25a) - (25b)                            //
       for (int t = 0; t < theta_.size(); 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]);
+        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]);
       }
     }
     double x2 = pow2(x.back());
@@ -1261,15 +1260,15 @@ c    MM       + 1  and - 1, alternately
       Qabs_ch_norm_[i] = Qext_ch_norm_[i] - Qsca_ch_norm_[i];
     }
     
-    albedo_ = Qsca_ / Qext_;                              // Equation (31)
-    asymmetry_factor_ = (Qext_ - Qpr_) / Qsca_;                          // Equation (32)
+    albedo_ = Qsca_/Qext_;                              // Equation (31)
+    asymmetry_factor_ = (Qext_ - Qpr_)/Qsca_;                          // Equation (32)
 
     Qbk_ = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2;    // Equation (33)
 
     isMieCalculated_ = true;
     nmax_used_ = nmax_;
     printf("Run Mie result: Qext = %g, Qsca = %g, Qabs = %g, Qbk = %g \n",
-    	   GetQext(), GetQsca(), GetQabs(), GetQbk() );
+               GetQext(), GetQsca(), GetQabs(), GetQbk());
     //return nmax;
   }
   
@@ -1284,10 +1283,10 @@ c    MM       + 1  and - 1, alternately
     //     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
-    al_n_.resize(L+1);
-    bl_n_.resize(L+1);
-    cl_n_.resize(L+1);
-    dl_n_.resize(L+1);
+    al_n_.resize(L + 1);
+    bl_n_.resize(L + 1);
+    cl_n_.resize(L + 1);
+    dl_n_.resize(L + 1);
     for (auto& element:al_n_) element.resize(nmax_);
     for (auto& element:bl_n_) element.resize(nmax_);
     for (auto& element:cl_n_) element.resize(nmax_);
@@ -1295,7 +1294,7 @@ c    MM       + 1  and - 1, alternately
     std::complex<double> c_one(1.0, 0.0);
     std::complex<double> c_zero(0.0, 0.0);
     // Yang, paragraph under eq. A3
-    // a^(L+1)_n = a_n, d^(L+1) = 1 ...
+    // a^(L + 1)_n = a_n, d^(L + 1) = 1 ...
     for (int i = 0; i < nmax_; ++i) {
       al_n_[L][i] = an_[i];
       bl_n_[L][i] = bn_[i];
@@ -1314,23 +1313,23 @@ c    MM       + 1  and - 1, alternately
     ScattCoeffsLayerdInit();
     const int L = index_.size();
     std::vector<std::complex<double> > z(L), z1(L);
-    for (int i = 0; i<L-1; ++i) {
+    for (int i = 0; i < L - 1; ++i) {
       z[i]  =size_parameter_[i]*index_[i];
-      z1[i]=size_parameter_[i]*index_[i+1];
+      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) {
-      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);
+      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);
     }
     for (int l = 0; l < L; ++l) {
       calcD1D3(z[l],D1z[l],D3z[l]);
@@ -1340,41 +1339,41 @@ c    MM       + 1  and - 1, alternately
     }
     auto& m = 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 (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 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]);
-	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]);
-	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]);
-	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]);
-	cl_n_[l][n] /= denom;   
+        // 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]);
+        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]);
+        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]);
+        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]);
+        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
@@ -1388,34 +1387,34 @@ c    MM       + 1  and - 1, alternately
     // 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());
+    //         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=%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) {
+    for (int i = 0; i < L + 1; ++i) {
       printf("Layer =%d ---> ", 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,
-    	       al_n_[i][n].real(), al_n_[i][n].imag(),
-	       bl_n_[i][n].real(), bl_n_[i][n].imag(),
-    	       cl_n_[i][n].real(), cl_n_[i][n].imag(),
-	       dl_n_[i][n].real(), dl_n_[i][n].imag());
+            //        if (n < 20) continue;
+            printf(" || n=%d --> a=%g,%g b=%g,%g c=%g,%g d=%g,%g",
+                   n,
+                   al_n_[i][n].real(), al_n_[i][n].imag(),
+                   bl_n_[i][n].real(), bl_n_[i][n].imag(),
+                   cl_n_[i][n].real(), cl_n_[i][n].imag(),
+                   dl_n_[i][n].real(), dl_n_[i][n].imag());
       }
       printf("\n\n");
     }
@@ -1431,15 +1430,14 @@ c    MM       + 1  and - 1, alternately
     
     std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0);
     std::vector<std::complex<double> > vm3o1n(3), vm3e1n(3), vn3o1n(3), vn3e1n(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> > bj(nmax_+1), by(nmax_+1), bd(nmax_+1);
+    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> > bj(nmax_ + 1), by(nmax_ + 1), bd(nmax_ + 1);
     // Calculate spherical Bessel and Hankel functions
     printf("##########  layer OUT ############\n");
     sphericalBessel(Rho,bj, by, bd);    
     for (int n = 0; n < nmax_; n++) {
       double rn = static_cast<double>(n + 1);
-      std::complex<double> zn = bj[n+1] + c_i*by[n+1];
+      std::complex<double> zn = bj[n + 1] + c_i*by[n + 1];
       // using BH 4.12 and 4.50
       std::complex<double> xxip = Rho*(bj[n] + c_i*by[n]) - rn*zn;
       
@@ -1461,10 +1459,10 @@ c    MM       + 1  and - 1, alternately
       // scattered field: BH p.94 (4.45)
       std::complex<double> encap = std::pow(c_i, rn)*(2.0*rn + 1.0)/(rn*rn + rn);
       for (int i = 0; i < 3; i++) {
-	Es[i] = Es[i] + encap*(c_i*an_[n]*vn3e1n[i] - bn_[n]*vm3o1n[i]);
-	Hs[i] = Hs[i] + encap*(c_i*bn_[n]*vn3o1n[i] + an_[n]*vm3e1n[i]);
-	//if (n<3) printf(" E[%d]=%g ", i,std::abs(Es[i]));
-	if (n<3) printf(" !!=%d=== %g ", i,std::abs(Es[i]));
+        Es[i] = Es[i] + encap*(c_i*an_[n]*vn3e1n[i] - bn_[n]*vm3o1n[i]);
+        Hs[i] = Hs[i] + encap*(c_i*bn_[n]*vn3o1n[i] + an_[n]*vm3e1n[i]);
+        //if (n<3) printf(" E[%d]=%g ", i,std::abs(Es[i]));
+        if (n<3) printf(" !!=%d=== %g ", i,std::abs(Es[i]));
       }
     }
     
@@ -1497,7 +1495,7 @@ c    MM       + 1  and - 1, alternately
     }
     
     for (int i = 0; i < 3; i++) {
-      // electric field E [V m-1] = EF*E0
+      // 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]));
@@ -1508,18 +1506,18 @@ c    MM       + 1  and - 1, alternately
   // ********************************************************************** //
   void MultiLayerMie::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)  {
     // printf("field int Qext = %g, Qsca = %g, Qabs = %g, Qbk = %g, \n",
-    // 	   GetQext(), GetQsca(), GetQabs(), GetQbk() );
+    //            GetQext(), GetQsca(), GetQabs(), GetQbk());
     
     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> > vm3o1n(3), vm3e1n(3), vn3o1n(3), vn3e1n(3);
     std::vector<std::complex<double> > vm1o1n(3), vm1e1n(3), vn1o1n(3), vn1e1n(3);
     std::vector<std::complex<double> > El(3,c_zero),Ei(3,c_zero), Hl(3,c_zero);
-    std::vector<std::complex<double> > bj(nmax_+1), by(nmax_+1), bd(nmax_+1);
+    std::vector<std::complex<double> > bj(nmax_ + 1), by(nmax_ + 1), bd(nmax_ + 1);
     int layer=0;  // layer number
     std::complex<double> index_l;
-    for (int i = 0; i < size_parameter_.size()-1; ++i) {
-      if (size_parameter_[i] < Rho && Rho <= size_parameter_[i+1]) {
-	layer=i;
+    for (int i = 0; i < size_parameter_.size() - 1; ++i) {
+      if (size_parameter_[i] < Rho && Rho <= size_parameter_[i + 1]) {
+        layer=i;
       }
     }
     if (Rho > size_parameter_.back()) {
@@ -1542,7 +1540,7 @@ c    MM       + 1  and - 1, alternately
     for (int n = 0; n < nmax_; n++) {
       double rn = static_cast<double>(n + 1);
       std::complex<double> znm1 = bj[n] + c_i*by[n];
-      std::complex<double> zn = bj[n+1] + c_i*by[n+1];
+      std::complex<double> zn = bj[n + 1] + c_i*by[n + 1];
       //if (n<3) printf("\nbesselh = %g,%g", zn.real(), zn.imag()); //!
       // using BH 4.12 and 4.50
       std::complex<double> xxip = Rho*(bj[n] + c_i*by[n]) - rn*zn;
@@ -1554,8 +1552,8 @@ c    MM       + 1  and - 1, alternately
       vm3o1n[1] = cos(Phi)*Pi[n]*zn;
       vm3o1n[2] = -sin(Phi)*Tau[n]*zn;
       // if (n<3)  printf("\nRE  vm3o1n[0]%g   vm3o1n[1]%g    vm3o1n[2]%g   \nIM vm3o1n[0]%g   vm3o1n[1]%g    vm3o1n[2]%g",
-      // 	     vm3o1n[0].real(),   vm3o1n[1].real(),    vm3o1n[2].real(),
-      // 	     vm3o1n[0].imag(),   vm3o1n[1].imag(),    vm3o1n[2].imag());
+      //              vm3o1n[0].real(),   vm3o1n[1].real(),    vm3o1n[2].real(),
+      //              vm3o1n[0].imag(),   vm3o1n[1].imag(),    vm3o1n[2].imag());
       vm3e1n[0] = c_zero;
       vm3e1n[1] = -sin(Phi)*Pi[n]*zn;
       vm3e1n[2] = -cos(Phi)*Tau[n]*zn;
@@ -1566,13 +1564,13 @@ c    MM       + 1  and - 1, alternately
       vn3e1n[1] = cos(Phi)*Tau[n]*xxip/Rho;
       vn3e1n[2] = -sin(Phi)*Pi[n]*xxip/Rho;
       // if (n<3)  printf("\nRE  vn3e1n[0]%g   vn3e1n[1]%g    vn3e1n[2]%g   \nIM vn3e1n[0]%g   vn3e1n[1]%g    vn3e1n[2]%g",
-      // 	     vn3e1n[0].real(),   vn3e1n[1].real(),    vn3e1n[2].real(),
-      // 	     vn3e1n[0].imag(),   vn3e1n[1].imag(),    vn3e1n[2].imag());
+      //              vn3e1n[0].real(),   vn3e1n[1].real(),    vn3e1n[2].real(),
+      //              vn3e1n[0].imag(),   vn3e1n[1].imag(),    vn3e1n[2].imag());
       
       znm1 = bj[n];
-      zn = bj[n+1];
+      zn = bj[n + 1];
       // znm1 = (bj[n] + c_i*by[n]).real();
-      // zn = (bj[n+1] + c_i*by[n+1]).real();
+      // zn = (bj[n + 1] + c_i*by[n + 1]).real();
       xxip = Rho*(bj[n]) - rn*zn;
       if (n<3)printf("\nbesselj = %g,%g", zn.real(), zn.imag()); //!
       vm1o1n[0] = c_zero;
@@ -1585,38 +1583,38 @@ c    MM       + 1  and - 1, alternately
       vn1o1n[1] = sin(Phi)*Tau[n]*xxip/Rho;
       vn1o1n[2] = cos(Phi)*Pi[n]*xxip/Rho;
       // if (n<3) printf("\nvn1o1n[2](%g) = cos(Phi)(%g)*Pi[n](%g)*xxip(%g)/Rho(%g)",
-      // 		      std::abs(vn1o1n[2]), cos(Phi),Pi[n],std::abs(xxip),Rho);
+      //                       std::abs(vn1o1n[2]), cos(Phi),Pi[n],std::abs(xxip),Rho);
       vn1e1n[0] = cos(Phi)*rn*(rn + 1.0)*sin(Theta)*Pi[n]*zn/Rho;
       vn1e1n[1] = cos(Phi)*Tau[n]*xxip/Rho;
       vn1e1n[2] = -sin(Phi)*Pi[n]*xxip/Rho;
       // if (n<3)  printf("\nRE  vm3o1n[0]%g   vm3o1n[1]%g    vm3o1n[2]%g   \nIM vm3o1n[0]%g   vm3o1n[1]%g    vm3o1n[2]%g",
-      // 	     vm3o1n[0].real(),   vm3o1n[1].real(),    vm3o1n[2].real(),
-      // 	     vm3o1n[0].imag(),   vm3o1n[1].imag(),    vm3o1n[2].imag());
+      //              vm3o1n[0].real(),   vm3o1n[1].real(),    vm3o1n[2].real(),
+      //              vm3o1n[0].imag(),   vm3o1n[1].imag(),    vm3o1n[2].imag());
       
       // scattered field: BH p.94 (4.45)
       std::complex<double> encap = std::pow(c_i, rn)*(2.0*rn + 1.0)/(rn*rn + rn);
       // if (n<3) printf("\n===== n=%d ======\n",n);
       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]
-		       );
-	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]
-			       + c_i*bl_n_[l][n]*vn3o1n[i] + al_n_[l][n]*vm3e1n[i]);
-	// printf("\n !Ei[%d]=%g,%g! ", i, Ei[i].real(), Ei[i].imag());
-	// if (n<3) printf("\n !El[%d]=%g,%g! ", i, El[i].real(), El[i].imag());
-	// //printf(" ===%d=== %g ", i,std::abs(cl_n_[l][n]*vm1o1n[i] - c_i*dl_n_[l][n]*vn1e1n[i]));
-	// if (n<3) printf(" ===%d=== %g ", i,std::abs(//-dl_n_[l][n]*vm1e1n[i] 
-	// 					    //- c_i*cl_n_[l][n]*
-	// 					    vn1o1n[i]
-	// 					    // + c_i*bl_n_[l][n]*vn3o1n[i]
-	// 					    // + al_n_[l][n]*vm3e1n[i]
-	// 					    ));
-	// if (n<3) printf(" --- Ei[%d]=%g! ", i,std::abs(encap*(vm1o1n[i] - c_i*vn1e1n[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]
+);
+        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]
+                               + c_i*bl_n_[l][n]*vn3o1n[i] + al_n_[l][n]*vm3e1n[i]);
+        // printf("\n !Ei[%d]=%g,%g! ", i, Ei[i].real(), Ei[i].imag());
+        // if (n<3) printf("\n !El[%d]=%g,%g! ", i, El[i].real(), El[i].imag());
+        // //printf(" ===%d=== %g ", i,std::abs(cl_n_[l][n]*vm1o1n[i] - c_i*dl_n_[l][n]*vn1e1n[i]));
+        // if (n<3) printf(" ===%d=== %g ", i,std::abs(//-dl_n_[l][n]*vm1e1n[i] 
+        //                                             //- c_i*cl_n_[l][n]*
+        //                                             vn1o1n[i]
+        //                                             // + c_i*bl_n_[l][n]*vn3o1n[i]
+        //                                             // + al_n_[l][n]*vm3e1n[i]
+        //                      ));
+        // if (n<3) printf(" --- Ei[%d]=%g! ", i,std::abs(encap*(vm1o1n[i] - c_i*vn1e1n[i])));
 
       }
       //if (n<3) printf(" bj=%g \n", std::abs(bj[n]));
@@ -1629,7 +1627,7 @@ c    MM       + 1  and - 1, alternately
     }
     
     for (int i = 0; i < 3; i++) {
-      // electric field E [V m-1] = EF*E0
+      // electric field E [V m - 1] = EF*E0
       E[i] = El[i];
       H[i] = Hl[i];
       printf("\n !El[%d]=%g,%g! ", i, El[i].real(), El[i].imag());
@@ -1647,8 +1645,8 @@ c    MM       + 1  and - 1, alternately
   // 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]     //
+  //   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.       //
@@ -1675,7 +1673,7 @@ c    MM       + 1  and - 1, alternately
     for (auto& f:E_field_) f.resize(3);
     for (auto& f:H_field_) f.resize(3);
 
-    for (int point = 0;	 point < total_points; ++point) {
+    for (int point = 0; point < total_points; ++point) {
       const double& Xp = coords_sp_[0][point];
       const double& Yp = coords_sp_[1][point];
       const double& Zp = coords_sp_[2][point];
@@ -1685,7 +1683,7 @@ c    MM       + 1  and - 1, alternately
       Rho = std::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
       // printf("Rho=%g\n", Rho);
       // Avoid convergence problems due to Rho too small
-      if (Rho < 1e-5) Rho = 1e-10;
+      if (Rho < 1e-10) Rho = 1e-10;
       // If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
       if (Rho == 0.0) Theta = 0.0;
       else Theta = std::acos(Zp/Rho);
@@ -1707,11 +1705,11 @@ c    MM       + 1  and - 1, alternately
       // Firstly the easiest case: the field outside the particle
       printf("rho=%g, outer=%g  ", Rho, outer_size);
       if (Rho >= outer_size) {
-      fieldExt(Rho, Phi, Theta, Pi, Tau, 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);
+        fieldExt(Rho, Phi, Theta, Pi, Tau, 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, Pi, Tau, 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);
+        fieldInt(Rho, Phi, Theta, Pi, Tau, 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);
       }
       std::complex<double>& Ex = E_field_[point][0];
       std::complex<double>& Ey = E_field_[point][1];
@@ -1721,18 +1719,18 @@ c    MM       + 1  and - 1, alternately
       std::complex<double>& Hz = H_field_[point][2];
       //Now, convert the fields back to cartesian coordinates
       {
-	using std::sin;
-	using std::cos;
-	Ex = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
-	Ey = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
-	Ez = cos(Theta)*Es[0] - sin(Theta)*Es[1];
+        using std::sin;
+        using std::cos;
+        Ex = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
+        Ey = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
+        Ez = cos(Theta)*Es[0] - sin(Theta)*Es[1];
       
-	Hx = sin(Theta)*cos(Phi)*Hs[0] + cos(Theta)*cos(Phi)*Hs[1] - sin(Phi)*Hs[2];
-	Hy = sin(Theta)*sin(Phi)*Hs[0] + cos(Theta)*sin(Phi)*Hs[1] + cos(Phi)*Hs[2];
-	Hz = cos(Theta)*Hs[0] - sin(Theta)*Hs[1];
+        Hx = sin(Theta)*cos(Phi)*Hs[0] + cos(Theta)*cos(Phi)*Hs[1] - sin(Phi)*Hs[2];
+        Hy = sin(Theta)*sin(Phi)*Hs[0] + cos(Theta)*sin(Phi)*Hs[1] + cos(Phi)*Hs[2];
+        Hz = 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);
+             Rho);
     }  // end of for all field coordinates
     
   }  //  end of   void MultiLayerMie::RunFieldCalculations()

nmie-wrapper.h

@@ -183,33 +183,33 @@ namespace nmie {
     void Nstop();
     void Nmax(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::vector<std::complex<double> >& jnp, std::vector<std::complex<double> >& h1n,
+	            std::vector<std::complex<double> >& h1np);
     void sphericalBessel(std::complex<double> z, std::vector<std::complex<double> >& bj,
-			 std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd);
+			             std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd);
     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> 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> 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);
+				                 double Pi, double Tau);
     std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,
-				 double Pi, double Tau);
+				                 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);
+		             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::complex<double> calcD1confra(int N, const std::complex<double> z);
     void calcD1D3(std::complex<double> z,
-		  std::vector<std::complex<double> >& D1,
-		  std::vector<std::complex<double> >& D3);
+		          std::vector<std::complex<double> >& D1,
+		          std::vector<std::complex<double> >& D3);
     void calcSinglePiTau(const double& costheta, std::vector<double>& Pi,
-			 std::vector<double>& Tau);
-    void calcAllPiTau( std::vector< std::vector<double> >& Pi,
-		    std::vector< std::vector<double> >& Tau);
+			             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();

nmie.cc

@@ -67,7 +67,7 @@ int Nstop(double xL) {
 //**********************************************************************************//
 int Nmax(int L, int fl, int pl,
          std::vector<double> x,
-		 std::vector<std::complex<double> > m) {
+         std::vector<std::complex<double> > m) {
   int i, result, ri, riM1;
   result = Nstop(x[L - 1]);
   for (i = fl; i < L; i++) {
@@ -251,8 +251,8 @@ void sphericalBessel(std::complex<double> z, int nmax, std::vector<std::complex<
 // 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)  {
+              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;
@@ -343,8 +343,8 @@ void fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double
 
 // 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> 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;
@@ -354,8 +354,8 @@ std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::com
 
 // 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> 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;
@@ -365,14 +365,14 @@ std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::com
 
 // Calculates S1 - equation (25a)
 std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
-		                     double Pi, double Tau) {
+                             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) {
+                             double Pi, double Tau) {
 
   return calc_S1(n, an, bn, Tau, Pi);
 }
@@ -391,10 +391,10 @@ std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double
 //   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) {
+                 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;
 
@@ -420,8 +420,8 @@ void calcPsiZeta(double x, int nmax,
 //   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) {
+              std::vector<std::complex<double> >& D1,
+              std::vector<std::complex<double> >& D3) {
 
   int n;
   std::complex<double> nz, PsiZeta;
@@ -496,7 +496,7 @@ void calcPiTau(int nmax, double Theta, std::vector<double>& Pi, std::vector<doub
 //   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) {
+                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 //
@@ -721,7 +721,7 @@ 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 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;

tests/python/field.py

@@ -1,25 +1,38 @@
 #!/usr/bin/env python
-
-# This test case calculates the differential scattering
-# cross section from a Luneburg lens, as described in:
-# B. R. Johnson, Applied Optics 35 (1996) 3286-3296.
-
-# The Luneburg lens is a sphere of radius a, with a
-# radially-varying index of refraction, given by:
-# m(r) = [2 - (r/a)**1]**(1/2)
-
-# For the calculations, the Luneburg lens was approximated
-# as a multilayered sphere with 500 equally spaced layers.
-# The refractive index of each layer is defined to be equal to
-# m(r) at the midpoint of the layer: ml = [2 - (xm/xL)**1]**(1/2),
-# with xm = (xl-1 + xl)/2, for l = 1,2,...,L. The size
-# parameter in the lth layer is xl = l*xL/500. According to
-# geometrical optics theory, the differential cross section
-# can be expressed as:
-# d(Csca)/d(a**2*Omega) = cos(Theta)
-
-# The differential cross section from wave optics is:
-# d(Csca)/d(a**2*Omega) = S11(Theta)/x**2
+# -*- 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 in the 
+# XY plane, 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

tests/python/lfield.py

@@ -0,0 +1,98 @@
+#!/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 the 
+# X, Y and Z axes, 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, 1), dtype = np.float64)
+x[0, 0] = 1.
+
+m = np.ones((1, 1), dtype = np.complex128)
+m[0, 0] = (0.05 + 2.070j)/1.46
+
+nc = 1001
+
+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)
+one = np.ones(nc, dtype = np.float64)
+
+coordX[:, 0] = scan
+coordY[:, 1] = scan
+coordZ[:, 2] = scan
+
+terms, Ex, Hx = fieldnlay(x, m, coordX)
+terms, Ey, Hy = fieldnlay(x, m, coordY)
+terms, Ez, Hz = fieldnlay(x, m, coordZ)
+
+Exr = np.absolute(Ex)
+Eyr = np.absolute(Ey)
+Ezr = np.absolute(Ez)
+
+# |E|/|Eo|
+Exh = np.sqrt(Exr[0, :, 0]**2 + Exr[0, :, 1]**2 + Exr[0, :, 2]**2)
+Eyh = np.sqrt(Eyr[0, :, 0]**2 + Eyr[0, :, 1]**2 + Eyr[0, :, 2]**2)
+Ezh = np.sqrt(Ezr[0, :, 0]**2 + Ezr[0, :, 1]**2 + Ezr[0, :, 2]**2)
+
+result = np.vstack((scan, Exh, Eyh, Ezh)).transpose()
+
+try:
+    import matplotlib.pyplot as plt
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+
+    ax.errorbar(result[:, 0], one, fmt = 'k')
+    ax.errorbar(result[:, 0], result[:, 1], fmt = 'r', label = 'X axis')
+    ax.errorbar(result[:, 0], result[:, 2], fmt = 'g', label = 'Y axis')
+    ax.errorbar(result[:, 0], result[:, 3], fmt = 'b', label = 'Z axis')
+
+    ax.legend()
+
+    plt.xlabel('X|Y|Z')
+    plt.ylabel('|E|/|Eo|')
+
+    plt.draw()
+    plt.show()
+
+finally:
+    np.savetxt("lfield.txt", result, fmt = "%.5f")
+    print result
+
+
+

tests/python/scattPEC.py

@@ -1,4 +1,29 @@
 #!/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 differential scattering
 # cross section for different x values of a PEC sphere

tests/python/test01.py

@@ -1,4 +1,29 @@
 #!/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 is a test against the program n-mie (version 3a) for the test case
 # distributed by them (extended for x up to 100)

tests/python/test02.py

@@ -1,4 +1,29 @@
 #!/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 is a test against the program n-mie (version 3a) for a water sphere surrounded by soot
 # n-mie is based in the algorithm described in:

tests/python/test04.py

@@ -1,4 +1,29 @@
 #!/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 differential scattering
 # cross section from a Luneburg lens, as described in: