stop treating access by reference as part of type (so changed 'int& i' to 'int &i')

examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -39,22 +39,25 @@ from matplotlib import pyplot as plt
 import inspect
 print("Using scattnlay from ", inspect.getfile(scattnlay))
 
-npts = 251
-factor = 3. # plot extent compared to sphere radius
-index_H2O = 1.33+0.j
+npts = 151
+factor = 1.50 # plot extent compared to sphere radius
+index_H2O = 1.33+(1e-6)*1j
 
 WL = 0.532 #mkm
-total_r = 1 #mkm
+total_r = 100 #mkm
 isMP = False
 # isMP = True
 
-# nmax = 230
-nmax = -1
 
 nm = 1.0  # host medium
 x = 2.0 * np.pi * np.array([total_r], dtype=np.float64) / WL
 m = np.array((index_H2O), dtype=np.complex128) / nm
 
+nmax = x*factor + 11 * (x*factor)**(1.0 / 3.0) + 1
+# return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
+
+# nmax = -1
+
 print("x =", x)
 print("m =", m)
 terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(

examples/example-eval-force.cc

@@ -92,7 +92,7 @@ int main(int argc, char *argv[]) {
       }
 
     }  // end for refines
-  } 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;
     return -1;

examples/example-get-Mie.cc

@@ -46,7 +46,7 @@ int main(int argc, char *argv[]) {
     //std::string TiN_filename("Si.txt");
     std::string shell_filename(core_filename);
 
-    nmie::MultiLayerMieApplied<nmie::FloatType> multi_layer_mie;  
+    nmie::MultiLayerMieApplied<nmie::FloatType> multi_layer_mie;
     const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
     const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
     const std::complex<double> index_Si = std::sqrt(epsilon_Si);
@@ -80,7 +80,7 @@ int main(int argc, char *argv[]) {
 
     multi_layer_mie.SetWavelength(WL);
     multi_layer_mie.RunMieCalculation();
-    
+
     double Qabs = static_cast<double>(multi_layer_mie.GetQabs());
     printf("Qabs = %g\n", Qabs);
     std::vector< std::vector<std::complex<nmie::FloatType> > > aln, bln, cln, dln;
@@ -110,11 +110,11 @@ int main(int argc, char *argv[]) {
     auto Si_data = Si_index.GetIndex();
     auto Ag_data = Ag_index.GetIndex();
     for (int i=0; i < Si_data.size(); ++i) {
-      const double& WL = Si_data[i].first;
-      const std::complex<double>& Si = Si_data[i].second;
-      const std::complex<double>& Ag = Ag_data[i].second;
+      const double &WL = Si_data[i].first;
+      const std::complex<double> &Si = Si_data[i].second;
+      const std::complex<double> &Ag = Ag_data[i].second;
       str+=std::to_string(WL);
-      multi_layer_mie.ClearTarget();      
+      multi_layer_mie.ClearTarget();
       if (isSiAgSi) {
 	multi_layer_mie.AddTargetLayer(core_width, Si);
 	multi_layer_mie.AddTargetLayer(inner_width, Ag);
@@ -202,7 +202,7 @@ int main(int argc, char *argv[]) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
-  }  
+  }
     return 0;
 }
 

examples/example-minimal.cc

@@ -52,7 +52,7 @@ int main(int , char **) {
     multi_layer_mie.RunMieCalculation();
     double Qabs = multi_layer_mie.GetQabs();
     printf("Qabs = %g\n", Qabs);
-  } 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;
     return -1;

examples/read-spectra.cc

@@ -2,11 +2,11 @@
  * @file   read-spectra.cc
  * @author Konstantin Ladutenko <kostyfisik at gmail (.) com>
  * @date   Wed Mar 11 11:51:26 2015
- * 
+ *
  * @copyright 2015 Konstantin Ladutenko
  *
  * @brief  Read complex spectra from file in format 'WL real imag'
- * 
+ *
  * read-spectra 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
@@ -19,7 +19,7 @@
  *
  * You should have received a copy of the GNU General Public License
  * along with read-spectra.  If not, see <http://www.gnu.org/licenses/>.
- * 
+ *
  */
 
 #include <algorithm>
@@ -37,7 +37,7 @@ namespace read_spectra {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  ReadSpectra& ReadSpectra::ReadFromFile(std::string fname) {
+  ReadSpectra &ReadSpectra::ReadFromFile(std::string fname) {
     //std::cout<<"Reading file: "<< fname << std::endl;
     std::ifstream infile(fname.c_str());
     data_.clear();
@@ -45,16 +45,16 @@ namespace read_spectra {
     while (std::getline(infile, line))
       {
 	if (line.front() == '#') continue; //do not read comments
-	if (line.find('#') != std::string::npos) 
+	if (line.find('#') != std::string::npos)
 	  throw std::invalid_argument("Error! Comments should be marked with # in the begining of the line!\n");
-	std::istringstream iss(line);	
+	std::istringstream iss(line);
 	double wl, re, im;
 	if (!(iss >> wl >> re >> im)) throw std::invalid_argument("Error! Unexpected format of the line!\n");
 	data_.push_back(std::vector<double>({wl,re,im}));
 	//std::cout<<wl<<' '<<re<<' '<<im<<std::endl;
-      }  // end of wile reading file 
+      }  // end of wile reading file
     std::sort(data_.begin(), data_.end(),
-	      [](const std::vector<double>& a, const std::vector<double>& b) {
+	      [](const std::vector<double> &a, const std::vector<double> &b) {
 		return a.front() < b.front();
 	      });
     return *this;
@@ -63,7 +63,7 @@ namespace read_spectra {
   // ********************************************************************** //
   // ********************************************************************** //
   /// Cut the spectra to the range and convert it to std::complex<double>
-  ReadSpectra& ReadSpectra::ResizeToComplex(double from_wl, double to_wl, int samples) {
+  ReadSpectra &ReadSpectra::ResizeToComplex(double from_wl, double to_wl, int samples) {
     if (data_.size() < 2) throw std::invalid_argument("Nothing to resize!/n");
     if (data_.front()[0] > from_wl || data_.front()[0] > to_wl ||
 	data_.back()[0] < from_wl || data_.back()[0] < to_wl ||
@@ -81,26 +81,26 @@ namespace read_spectra {
     data_complex_.clear();
     int j = 0;
     for (int i = 0; i < data_.size(); ++i) {
-      const double& wl_i = data_[i][0];
-      const double& wl_s = wl_sampled[j];
+      const double &wl_i = data_[i][0];
+      const double &wl_s = wl_sampled[j];
       if (wl_i < wl_s) continue;
       else {
-	const double& wl_prev = data_[i-1][0];	
-	const double& re_prev = data_[i-1][1];
-	const double& im_prev = data_[i-1][2];
-	const double& re_i = data_[i][1];
-	const double& im_i = data_[i][2];
+	const double &wl_prev = data_[i-1][0];
+	const double &re_prev = data_[i-1][1];
+	const double &im_prev = data_[i-1][2];
+	const double &re_i = data_[i][1];
+	const double &im_i = data_[i][2];
 	// Linear approximation
 	double re_s = re_i - (re_i-re_prev)*(wl_i-wl_s)/(wl_i-wl_prev);
 	double im_s = im_i - (im_i-im_prev)*(wl_i-wl_s)/(wl_i-wl_prev);
 
 	auto tmp = std::make_pair(wl_s, std::complex<double>(re_s,im_s));
-	data_complex_.push_back(tmp);	
-	       
+	data_complex_.push_back(tmp);
+
 	++j;
 	--i; // Next sampled point(j) can be in the same i .. i-1 region
 	// All sampled wavelengths has a value
-	if (j >= wl_sampled.size()) break;  
+	if (j >= wl_sampled.size()) break;
       }
     }
     if (data_complex_.size() == 0)
@@ -113,7 +113,7 @@ namespace read_spectra {
   // ********************************************************************** //
   // ********************************************************************** //
   /// from relative permittivity to refractive index
-  ReadSpectra& ReadSpectra::ToIndex() {
+  ReadSpectra &ReadSpectra::ToIndex() {
     data_complex_index_.clear();
     for (auto row : data_complex_) {
       const double wl = row.first;
@@ -122,7 +122,7 @@ namespace read_spectra {
       const double n = std::sqrt( (std::sqrt(pow2(e1)+pow2(e2)) + e1) /2.0 );
       const double k = std::sqrt( (std::sqrt(pow2(e1)+pow2(e2)) - e1) /2.0 );
       auto tmp = std::make_pair(wl, std::complex<double>(n,k));
-      data_complex_index_.push_back(tmp);	
+      data_complex_index_.push_back(tmp);
     }
     return *this;
   }
@@ -131,7 +131,7 @@ namespace read_spectra {
   // ********************************************************************** //
   // ********************************************************************** //
   void ReadSpectra::PrintData() {
-    if (data_complex_.size() == 0) 
+    if (data_complex_.size() == 0)
       throw std::invalid_argument("Nothing to print!");
     for (auto row : data_complex_) {
       printf("wl:%g\tre:%g\tim:%g\n", row.first, row.second.real(),

examples/read-spectra.h

@@ -7,7 +7,7 @@
  * @copyright 2015 Konstantin Ladutenko
  *
  * @brief  Read complex spectra from file in format 'WL real imag'
- * 
+ *
  * read-spectra 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
@@ -20,7 +20,7 @@
  *
  * You should have received a copy of the GNU General Public License
  * along with read-spectra.  If not, see <http://www.gnu.org/licenses/>.
- * 
+ *
  */
 #include <vector>
 #include <string>
@@ -29,11 +29,11 @@
 namespace read_spectra {
   class ReadSpectra {  // will throw for any error
    public:
-    ReadSpectra& ReadFromFile(std::string filename);
-    ReadSpectra& ResizeToComplex(double from_wl, double to_wl, int samples);
-    ReadSpectra& ToIndex();
+    ReadSpectra &ReadFromFile(std::string filename);
+    ReadSpectra &ResizeToComplex(double from_wl, double to_wl, int samples);
+    ReadSpectra &ToIndex();
     std::complex<double> at(double wavelength);
-    void PrintData();    
+    void PrintData();
     std::vector< std::pair< double, std::complex<double> > >&
       GetIndex(){return data_complex_index_;};
   private:

examples/test-surf-integral.cc

@@ -44,12 +44,12 @@ int main(int argc, char *argv[]) {
     //shell.PrintVerts();
     auto points = shell.GetVertices();
     double charge_s = shell.IntegrateGaussSimple(charge, shift);
-    std::cout << "Accuracy ( 1==ideal ): " << charge_s/charge << std::endl; 
-  } catch( const std::invalid_argument& ia ) {
+    std::cout << "Accuracy ( 1==ideal ): " << charge_s/charge << std::endl;
+  } catch( const std::invalid_argument &ia ) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;
     return -1;
-  }  
+  }
     return 0;
 }
 

src/farfield.cc

@@ -80,7 +80,7 @@ int main(int argc, char *argv[]) {
 
     int mode = -1;
     double tmp_mr;
-    for (const auto& arg : args) {
+    for (const auto &arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.

src/nearfield.cc

@@ -80,7 +80,7 @@ int main(int argc, char *argv[]) {
 
     int mode = -1;
     double tmp_mr;
-    for (const auto& arg : args) {
+    for (const auto &arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
@@ -252,7 +252,7 @@ int main(int argc, char *argv[]) {
     }
 
 
-  } 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;
     return -1;

src/nmie-applied.hpp

@@ -43,7 +43,7 @@
 
 namespace nmie {
 
-  int nMieApplied(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+  int nMieApplied(const unsigned int L, const int pl, std::vector<double> &x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMieApplied(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMieApplied(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMieApplied(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);

src/nmie-basic.hpp

@@ -163,7 +163,7 @@ namespace nmie {
   // Modify scattering (theta) angles                                       //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::SetAngles(const std::vector<FloatType>& angles) {
+  void MultiLayerMie<FloatType>::SetAngles(const std::vector<FloatType> &angles) {
     MarkUncalculated();
     theta_ = angles;
   }
@@ -173,7 +173,7 @@ namespace nmie {
   // Modify size of all layers                                             //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::SetLayersSize(const std::vector<FloatType>& layer_size) {
+  void MultiLayerMie<FloatType>::SetLayersSize(const std::vector<FloatType> &layer_size) {
     MarkUncalculated();
     size_param_.clear();
     FloatType prev_layer_size = 0.0;
@@ -193,7 +193,7 @@ namespace nmie {
   // Modify refractive index of all layers                                  //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::SetLayersIndex(const std::vector< std::complex<FloatType> >& index) {
+  void MultiLayerMie<FloatType>::SetLayersIndex(const std::vector< std::complex<FloatType> > &index) {
     MarkUncalculated();
     refractive_index_ = index;
   }
@@ -203,7 +203,7 @@ namespace nmie {
   // Modify coordinates for field calculation                               //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::SetFieldCoords(const std::vector< std::vector<FloatType> >& coords) {
+  void MultiLayerMie<FloatType>::SetFieldCoords(const std::vector< std::vector<FloatType> > &coords) {
     if (coords.size() != 3)
       throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
     if (coords[0].size() != coords[1].size() || coords[0].size() != coords[2].size())
@@ -312,8 +312,8 @@ namespace nmie {
     const int pl = PEC_layer_position_;
     const unsigned int first_layer = (pl > 0) ? pl : 0;
     unsigned int ri, riM1, nmax = 0;
-    const std::vector<FloatType>& x = size_param_;
-    const std::vector<std::complex<FloatType> >& m = refractive_index_;
+    const std::vector<FloatType> &x = size_param_;
+    const std::vector<std::complex<FloatType> > &m = refractive_index_;
     nmax = calcNstop(xL);
     for (unsigned int i = first_layer; i < x.size(); i++) {
       if (static_cast<int>(i) > PEC_layer_position_)  // static_cast used to avoid warning
@@ -408,8 +408,8 @@ namespace nmie {
   //**********************************************************************************//
   template <typename FloatType>
   void MultiLayerMie<FloatType>::calcD1D3(const std::complex<FloatType> z,
-                               std::vector<std::complex<FloatType> >& D1,
-                               std::vector<std::complex<FloatType> >& D3) {
+                               std::vector<std::complex<FloatType> > &D1,
+                               std::vector<std::complex<FloatType> > &D3) {
     std::vector<std::complex<FloatType> > PsiZeta(nmax_+1);
     evalDownwardD1(z, D1);
     evalUpwardD3 (z, D1, D3, PsiZeta);
@@ -430,8 +430,8 @@ namespace nmie {
   //**********************************************************************************//
   template <typename FloatType>
   void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
-                                  std::vector<std::complex<FloatType> >& Psi,
-                                  std::vector<std::complex<FloatType> >& Zeta) {
+                                  std::vector<std::complex<FloatType> > &Psi,
+                                  std::vector<std::complex<FloatType> > &Zeta) {
     std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1),
         PsiZeta(nmax_+1);
     // First, calculate the logarithmic derivatives
@@ -447,6 +447,22 @@ namespace nmie {
   }
 
 
+  template <typename FloatType>
+  void MultiLayerMie<FloatType>::calcPiTauAllTheta(const double from_Theta, const double to_Theta,
+                                                   std::vector<std::vector<FloatType> > &Pi,
+                                                   std::vector<std::vector<FloatType> > &Tau) {
+    auto perimeter_points = Pi.size();
+    for (auto &val:Pi) val.resize(available_maximal_nmax_);
+    for (auto &val:Tau) val.resize(available_maximal_nmax_);
+    double delta_Theta = eval_delta<FloatType>(perimeter_points, from_Theta, to_Theta);
+    for (unsigned int i=0; i < perimeter_points; i++) {
+      auto Theta = static_cast<FloatType>(from_Theta + i*delta_Theta);
+      // Calculate angular functions Pi and Tau
+      calcPiTau(nmm::cos(Theta), Pi[i], Tau[i]);
+    }
+  }
+
+
   //**********************************************************************************//
   // This function calculates Pi and Tau for a given value of cos(Theta).             //
   // Equations (26a) - (26c)                                                          //
@@ -461,8 +477,8 @@ namespace nmie {
   //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
   //**********************************************************************************//
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcPiTau(const FloatType& costheta,
-                                std::vector<FloatType>& Pi, std::vector<FloatType>& Tau) {
+  void MultiLayerMie<FloatType>::calcPiTau(const FloatType &costheta,
+                                std::vector<FloatType> &Pi, std::vector<FloatType> &Tau) {
 
     int nmax = Pi.size();
     if (Pi.size() != Tau.size())
@@ -504,10 +520,10 @@ namespace nmie {
   //**********************************************************************************//
   template <typename FloatType>
   void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
-                                    const std::complex<FloatType>& rn, const std::complex<FloatType>& Dn,
-                                    const FloatType& Pi, const FloatType& Tau, const FloatType& n,
-                                    std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n,
-                                    std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n) {
+                                    const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
+                                    const FloatType &Pi, const FloatType &Tau, const FloatType &n,
+                                    std::vector<std::complex<FloatType> > &Mo1n, std::vector<std::complex<FloatType> > &Me1n,
+                                    std::vector<std::complex<FloatType> > &No1n, std::vector<std::complex<FloatType> > &Ne1n) {
 
     // using eq 4.50 in BH
     std::complex<FloatType> c_zero(0.0, 0.0);
@@ -556,9 +572,9 @@ namespace nmie {
 
     isScaCoeffsCalc_ = false;
 
-    const std::vector<FloatType>& x = size_param_;
-    const std::vector<std::complex<FloatType> >& m = refractive_index_;
-    const int& pl = PEC_layer_position_;
+    const std::vector<FloatType> &x = size_param_;
+    const std::vector<std::complex<FloatType> > &m = refractive_index_;
+    const int &pl = PEC_layer_position_;
     const int L = refractive_index_.size();
 
 
@@ -753,7 +769,7 @@ namespace nmie {
     if (size_param_.size() == 0)
       throw std::invalid_argument("Initialize model first!");
 
-    const std::vector<FloatType>& x = size_param_;
+    const std::vector<FloatType> &x = size_param_;
 
     //MarkUncalculated();
 

src/nmie-js-wrapper.cc

@@ -85,7 +85,7 @@ EMSCRIPTEN_BINDINGS (c) {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-//  int nMieApplied(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+//  int nMieApplied(const unsigned int L, const int pl, std::vector<double> &x, std::vector<std::complex<double> > &m, const unsigned int nTheta, std::vector<double> &Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
 //
 //    if (x.size() != L || m.size() != L)
 //        throw std::invalid_argument("Declared number of layers do not fit x and m!");

src/nmie-nearfield.hpp

@@ -111,7 +111,7 @@ namespace nmie {
     std::vector<std::complex<FloatType> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
     std::complex<FloatType> denomZeta, denomPsi, T1, T2, T3, T4;
 
-    auto& m = refractive_index_;
+    auto &m = refractive_index_;
     std::vector< std::complex<FloatType> > m1(L);
 
     for (int l = 0; l < L - 1; l++) m1[l] = m[l + 1];
@@ -211,8 +211,8 @@ namespace nmie {
   template <typename FloatType>
   void MultiLayerMie<FloatType>::calcFieldByComponents(const FloatType Rho,
                                 const FloatType Theta, const FloatType Phi,
-                                std::vector<std::complex<FloatType> >& E,
-                                std::vector<std::complex<FloatType> >& H)  {
+                                std::vector<std::complex<FloatType> > &E,
+                                std::vector<std::complex<FloatType> > &H)  {
 
     std::complex<FloatType> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
     // Vector containing precomputed integer powers of i to avoid computation
@@ -222,7 +222,6 @@ namespace nmie {
     std::vector<std::complex<FloatType> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
     std::vector<FloatType> Pi(nmax_), Tau(nmax_);
 
-    int l = 0;  // Layer number
     std::complex<FloatType> ml;
 
     // Initialize E and H
@@ -231,17 +230,8 @@ namespace nmie {
       H[i] = c_zero;
     }
 
-    if (Rho > size_param_.back()) {
-      l = size_param_.size();
-      ml = c_one;
-    } else {
-      for (int i = size_param_.size() - 1; i >= 0 ; i--) {
-        if (Rho <= size_param_[i]) {
-          l = i;
-        }
-      }
-      ml = refractive_index_[l];
-    }
+    unsigned int l;
+    GetIndexAtRadius(Rho, ml, l);
 
     // Calculate logarithmic derivative of the Ricatti-Bessel functions
     calcD1D3(Rho*ml, D1n, D3n);
@@ -353,15 +343,15 @@ namespace nmie {
     H_.resize(total_points);
     Es_.resize(total_points);
     Hs_.resize(total_points);
-    for (auto& f : E_) f.resize(3);
-    for (auto& f : H_) f.resize(3);
-    for (auto& f : Es_) f.resize(3);
-    for (auto& f : Hs_) f.resize(3);
+    for (auto &f : E_) f.resize(3);
+    for (auto &f : H_) f.resize(3);
+    for (auto &f : Es_) f.resize(3);
+    for (auto &f : Hs_) f.resize(3);
 
     for (int point = 0; point < total_points; point++) {
-      const FloatType& Xp = coords_[0][point];
-      const FloatType& Yp = coords_[1][point];
-      const FloatType& Zp = coords_[2][point];
+      const FloatType &Xp = coords_[0][point];
+      const FloatType &Yp = coords_[1][point];
+      const FloatType &Zp = coords_[2][point];
 
       // Convert to spherical coordinates
       Rho = nmm::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
@@ -416,7 +406,9 @@ int ceil_to_2_pow_n(const int input_n) {
 
 template <typename FloatType>
 double eval_delta(const int steps, const double from_value, const double to_value) {
-  auto delta = std::abs((from_value - to_value)/static_cast<double>(steps));
+  auto delta = std::abs(from_value - to_value);
+  if (steps < 2) return delta;
+  delta /= static_cast<double>(steps-1);
   // We have a limited double precision evaluation of special functions, typically it is 1e-10.
   if ( (2.*delta)/std::abs(from_value+to_value) < 1e-9)
     throw std::invalid_argument("Error! The step is too fine, not supported!");
@@ -424,6 +416,78 @@ double eval_delta(const int steps, const double from_value, const double to_valu
 }
 
 
+// ml - refractive index
+// l - Layer number
+template <typename FloatType>
+void MultiLayerMie<FloatType>::GetIndexAtRadius(const FloatType Rho,
+                                                std::complex<FloatType> &ml,
+                                                unsigned int &l) {
+  l = 0;
+  if (Rho > size_param_.back()) {
+    l = size_param_.size();
+    ml = std::complex<FloatType>(1.0, 0.0);
+  } else {
+    for (int i = size_param_.size() - 1; i >= 0 ; i--) {
+      if (Rho <= size_param_[i]) {
+        l = i;
+      }
+    }
+    ml = refractive_index_[l];
+  }
+}
+template <typename FloatType>
+void MultiLayerMie<FloatType>::GetIndexAtRadius(const FloatType Rho,
+                                                std::complex<FloatType> &ml) {
+  unsigned int l;
+  GetIndexAtRadius(Rho, ml, l);
+}
+
+template <typename FloatType>
+void MultiLayerMie<FloatType>::calcMieSeriesNeededToConverge(const FloatType Rho) {
+  auto required_near_field_nmax = calcNmax(Rho);
+  SetMaxTerms(required_near_field_nmax);
+  // Calculate scattering coefficients an_ and bn_
+  calcScattCoeffs();
+  // We might be limited with available machine precision
+  available_maximal_nmax_ = nmax_;
+  // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
+  calcExpanCoeffs();
+}
+
+
+template <typename FloatType>
+void MultiLayerMie<FloatType>::calcRadialOnlyDependantFunctions(const FloatType from_Rho, const FloatType to_Rho,
+                                                                const bool isIgnoreAvailableNmax,
+                                                                std::vector<std::vector<std::complex<FloatType> > > &Psi,
+                                                                std::vector<std::vector<std::complex<FloatType> > > &D1n,
+                                                                std::vector<std::vector<std::complex<FloatType> > > &Zeta,
+                                                                std::vector<std::vector<std::complex<FloatType> > > &D3n) {
+  unsigned int radius_points = Psi.size();
+  std::vector<std::vector<std::complex<FloatType> > > PsiZeta(radius_points);
+  double delta_Rho = eval_delta<FloatType>(radius_points, from_Rho, to_Rho);
+  for (unsigned int j=0; j < radius_points; j++) {
+    auto Rho = static_cast<FloatType>(from_Rho + j*delta_Rho);
+    int near_field_nmax = calcNmax(to_Rho);
+    // Skip if not enough terms in Mie series (i.e. required near field nmax > available terms )
+    if (near_field_nmax > available_maximal_nmax_ && !isIgnoreAvailableNmax) continue;
+    if (near_field_nmax > available_maximal_nmax_)  near_field_nmax = available_maximal_nmax_;
+    Psi[j].resize(near_field_nmax + 1); D1n[j].resize(near_field_nmax + 1);
+    Zeta[j].resize(near_field_nmax + 1); D3n[j].resize(near_field_nmax + 1);
+    PsiZeta[j].resize(near_field_nmax + 1);
+    std::complex<FloatType> ml;
+    GetIndexAtRadius(Rho, ml);
+    auto z = Rho*ml;
+    evalDownwardD1(z, D1n[j]);
+    evalUpwardPsi(z,  D1n[j], Psi[j]);
+    evalUpwardD3 (z, D1n[j], D3n[j], PsiZeta[j]);
+    for (unsigned int k = 0; k < Zeta[j].size(); k++) {
+      Zeta[j][k] = PsiZeta[j][k]/Psi[j][k];
+    }
+  }
+
+}
+
+
 // input parameters:
 //         input_outer_perimeter_points: will be increased to the nearest power of 2.
 template <typename FloatType>
@@ -431,33 +495,97 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int input_outer_pe
                                                         const int radius_points,
                                                         const double from_Rho, const double to_Rho,
                                                         const double from_Theta, const double to_Theta,
-                                                        const double from_Phi, const double to_Phi) {
+                                                        const double from_Phi, const double to_Phi,
+                                                        const bool isIgnoreAvailableNmax) {
 //  double Rho, Theta, Phi;
   if (from_Rho > to_Rho || from_Theta > to_Theta || from_Phi > to_Phi
-      || input_outer_perimeter_points < 1 || radius_points < 1)
+      || input_outer_perimeter_points < 1 || radius_points < 1
+      || from_Rho < 0.)
     throw std::invalid_argument("Error! Invalid argument for RunFieldCalculationPolar() !");
   int outer_perimeter_points = input_outer_perimeter_points;
   if (outer_perimeter_points != 1) outer_perimeter_points = ceil_to_2_pow_n<FloatType>(input_outer_perimeter_points);
-//  double delta_Rho = eval_delta<FloatType>(radius_points, from_Rho, to_Rho);
-//  double delta_Phi = eval_delta<FloatType>(radius_points, from_Phi, to_Phi);
-  double delta_Theta = eval_delta<FloatType>(radius_points, from_Theta, to_Theta);
-  auto near_field_nmax = calcNmax(to_Rho);
-  SetMaxTerms(near_field_nmax);
+
+  calcMieSeriesNeededToConverge(to_Rho);
 
   std::vector<std::vector<FloatType> >  Pi(outer_perimeter_points), Tau(outer_perimeter_points);
-  for (auto &val:Pi) val.resize(near_field_nmax);
-  for (auto &val:Tau) val.resize(near_field_nmax);
-  for (int i=0; i < outer_perimeter_points; i++) {
-    auto Theta = static_cast<FloatType>(from_Theta + i*delta_Theta);
-    // Calculate angular functions Pi and Tau
-    calcPiTau(nmm::cos(Theta), Pi[i], Tau[i]);
-  }
+  calcPiTauAllTheta(from_Theta, to_Theta, Pi, Tau);
 
+  std::vector<std::vector<std::complex<FloatType> > > Psi(radius_points), D1n(radius_points),
+      Zeta(radius_points), D3n(radius_points), PsiZeta(radius_points);
+  calcRadialOnlyDependantFunctions(from_Rho, to_Rho, isIgnoreAvailableNmax,
+                                   Psi, D1n, Zeta, D3n);
+
+//  double delta_Phi = eval_delta<FloatType>(radius_points, from_Phi, to_Phi);
+  //  std::complex<FloatType> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
+//  // Vector containing precomputed integer powers of i to avoid computation
+//  std::vector<std::complex<FloatType> > ipow = {c_one, c_i, -c_one, -c_i};
+//  std::vector<std::complex<FloatType> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
+//  std::vector<std::complex<FloatType> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
+//  std::vector<std::complex<FloatType> > E, H
+//  std::complex<FloatType> ml;
+
+//  // Initialize E and H
+//  for (int i = 0; i < 3; i++) {
+//    E[i] = c_zero;
+//    H[i] = c_zero;
+//  }
+//
+//  auto ml = GetIndexAtRadius(Rho);
+//
+//  for (int n = 0; n < nmax_; n++) {
+//    int n1 = n + 1;
+//    auto rn = static_cast<FloatType>(n1);
+//
+//    // using BH 4.12 and 4.50
+//    calcSpherHarm(Rho*ml, Theta, Phi, Psi[n1], D1n[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
+//    calcSpherHarm(Rho*ml, Theta, Phi, Zeta[n1], D3n[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<FloatType> En = ipow[n1 % 4]
+//        *static_cast<FloatType>((rn + rn + 1.0)/(rn*rn + rn));
+//    for (int i = 0; i < 3; i++) {
+//      auto Ediff = 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]);
+//      auto Hdiff = 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]);
+//      if (nmm::isnan(Ediff.real()) || nmm::isnan(Ediff.imag()) ||
+//          nmm::isnan(Hdiff.real()) || nmm::isnan(Hdiff.imag())
+//          ) break;
+//      if (mode_n_ == Modes::kAll) {
+//        // electric field E [V m - 1] = EF*E0
+//        E[i] += Ediff;
+//        H[i] += Hdiff;
+//        continue;
+//      }
+//      if (n1 == mode_n_) {
+//        if (mode_type_ == Modes::kElectric || mode_type_ == Modes::kAll) {
+//          E[i] += En*( -c_i*dln_[l][n]*N1e1n[i]
+//              + c_i*aln_[l][n]*N3e1n[i]);
+//
+//          H[i] += En*(-dln_[l][n]*M1e1n[i]
+//              +aln_[l][n]*M3e1n[i]);
+//          //std::cout << mode_n_;
+//        }
+//        if (mode_type_ == Modes::kMagnetic  || mode_type_ == Modes::kAll) {
+//          E[i] += En*(  cln_[l][n]*M1o1n[i]
+//              - bln_[l][n]*M3o1n[i]);
+//
+//          H[i] += En*( -c_i*cln_[l][n]*N1o1n[i]
+//              + c_i*bln_[l][n]*N3o1n[i]);
+//          //std::cout << mode_n_;
+//        }
+//        //std::cout << std::endl;
+//      }
+//      //throw std::invalid_argument("Error! Unexpected mode for field evaluation!\n mode_n="+std::to_string(mode_n)+", mode_type="+std::to_string(mode_type)+"\n=====*****=====");
+//    }
+//  }  // end of for all n
+//
+//  // magnetic field
+//  std::complex<FloatType> hffact = ml/static_cast<FloatType>(cc_*mu_);
+//  for (int i = 0; i < 3; i++) {
+//    H[i] = hffact*H[i];
+//  }
 
-  // Calculate scattering coefficients an_ and bn_
-//  calcScattCoeffs();
-  // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
-//  calcExpanCoeffs();
 }
 }  // end of namespace nmie
 #endif  // SRC_NMIE_NEARFIELD_HPP_

src/nmie-pybind11.cc

@@ -50,13 +50,13 @@ py::array_t< std::complex<double>> VectorComplex2Py(const std::vector<std::compl
 
 // https://stackoverflow.com/questions/17294629/merging-flattening-sub-vectors-into-a-single-vector-c-converting-2d-to-1d
 template <typename T>
-std::vector<T> flatten(const std::vector<std::vector<T>>& v) {
+std::vector<T> flatten(const std::vector<std::vector<T>> &v) {
     std::size_t total_size = 0;
-    for (const auto& sub : v)
+    for (const auto &sub : v)
         total_size += sub.size(); // I wish there was a transform_accumulate
     std::vector<T> result;
     result.reserve(total_size);
-    for (const auto& sub : v)
+    for (const auto &sub : v)
         result.insert(result.end(), sub.begin(), sub.end());
     return result;
 }
@@ -162,8 +162,8 @@ py::tuple py_fieldnlay(const py::array_t<double, py::array::c_style | py::array:
   unsigned int ncoord = py_Xp.size();
   std::vector<std::vector<std::complex<double> > > E(ncoord);
   std::vector<std::vector<std::complex<double> > > H(ncoord);
-  for (auto& f : E) f.resize(3);
-  for (auto& f : H) f.resize(3);
+  for (auto &f : E) f.resize(3);
+  for (auto &f : H) f.resize(3);
   int L = py_x.size(), terms;
   terms = nmie::nField(L, pl, c_x, c_m, nmax, nmie::Modes::kAll, nmie::Modes::kAll, ncoord, c_Xp, c_Yp, c_Zp, E, H);
   auto py_E = Vector2DComplex2Py<std::complex<double> >(E);

src/nmie.cc

@@ -74,8 +74,8 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ScattCoeffs(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m,
-                  const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
+  int ScattCoeffs(const unsigned int L, const int pl, const std::vector<double> &x, const std::vector<std::complex<double> > &m,
+                  const int nmax, std::vector<std::complex<double> > &an, std::vector<std::complex<double> > &bn) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -92,7 +92,7 @@ namespace nmie {
       bn = ConvertComplexVector<double>(ml_mie.GetBn());
 
       return ml_mie.GetMaxTerms();
-    } catch(const std::invalid_argument& ia) {
+    } catch(const std::invalid_argument &ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
@@ -118,9 +118,9 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ExpanCoeffs(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax,
-                  std::vector<std::vector<std::complex<double> > >& aln, std::vector<std::vector<std::complex<double> > >& bln,
-                  std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
+  int ExpanCoeffs(const unsigned int L, const int pl, const std::vector<double> &x, const std::vector<std::complex<double> > &m, const int nmax,
+                  std::vector<std::vector<std::complex<double> > > &aln, std::vector<std::vector<std::complex<double> > > &bln,
+                  std::vector<std::vector<std::complex<double> > > &cln, std::vector<std::vector<std::complex<double> > > &dln) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -142,7 +142,7 @@ namespace nmie {
       dln = ConvertComplexVectorVector<double>(ml_mie.GetLayerDn());
 
       return ml_mie.GetMaxTerms();
-    } catch(const std::invalid_argument& ia) {
+    } catch(const std::invalid_argument &ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
@@ -179,10 +179,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta, const int nmax,
+  int nMie(const unsigned int L, const int pl, std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta, const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2,
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2,
            int mode_n, int mode_type) {
 
     if (x.size() != L || m.size() != L)
@@ -217,7 +217,7 @@ namespace nmie {
       S2 = ConvertComplexVector<double>(ml_mie.GetS2());
 
       return ml_mie.GetMaxTerms();
-    } catch(const std::invalid_argument& ia) {
+    } catch(const std::invalid_argument &ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
@@ -256,10 +256,10 @@ namespace nmie {
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
   int nMie(const unsigned int L, const int pl,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta, const int nmax,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta, const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2) {
     return nmie::nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, -1, -1);
   }
   //**********************************************************************************//
@@ -289,10 +289,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+  int nMie(const unsigned int L, std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2) {
     return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -324,10 +324,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+  int nMie(const unsigned int L, const int pl, std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2) {
     return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -361,10 +361,10 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta, const int nmax,
+  int nMie(const unsigned int L, std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta, const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2) {
     return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -392,19 +392,19 @@ namespace nmie {
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
   int nField(const unsigned int L, const int pl,
-             const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax,
+             const std::vector<double> &x, const std::vector<std::complex<double> > &m, const int nmax,
              const int mode_n, const int mode_type, const unsigned 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) {
+             const std::vector<double> &Xp_vec, const std::vector<double> &Yp_vec, const std::vector<double> &Zp_vec,
+             std::vector<std::vector<std::complex<double> > > &E, std::vector<std::vector<std::complex<double> > > &H) {
     if (x.size() != L || m.size() != L)
       throw std::invalid_argument("Declared number of layers do not fit x and m!");
     if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
         || E.size() != ncoord || H.size() != ncoord)
       throw std::invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
-    for (const auto& f:E)
+    for (const auto &f:E)
       if (f.size() != 3)
         throw std::invalid_argument("Field E is not 3D!");
-    for (const auto& f:H)
+    for (const auto &f:H)
       if (f.size() != 3)
         throw std::invalid_argument("Field H is not 3D!");
     try {
@@ -423,7 +423,7 @@ namespace nmie {
       H = ConvertComplexVectorVector<double>(ml_mie.GetFieldH());
 
       return ml_mie.GetMaxTerms();
-    } catch(const std::invalid_argument& ia) {
+    } catch(const std::invalid_argument &ia) {
       // Will catch if  ml_mie fails or other errors.
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);

src/nmie.hpp

@@ -43,20 +43,22 @@
 #endif
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl,
-                  const std::vector<double>& x, const std::vector<std::complex<double> >& m,
+                  const std::vector<double> &x, const std::vector<std::complex<double> > &m,
                   const 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);
 
   int ExpanCoeffs(const unsigned int L, const int pl,
-                  const std::vector<double>& x, const std::vector<std::complex<double> >& m,
+                  const std::vector<double> &x, const std::vector<std::complex<double> > &m,
                   const int nmax,
-                  std::vector<std::vector<std::complex<double> > >& an,
-                  std::vector<std::vector<std::complex<double> > >& bn,
-                  std::vector<std::vector<std::complex<double> > >& cn,
-                  std::vector<std::vector<std::complex<double> > >& dn);
+                  std::vector<std::vector<std::complex<double> > > &an,
+                  std::vector<std::vector<std::complex<double> > > &bn,
+                  std::vector<std::vector<std::complex<double> > > &cn,
+                  std::vector<std::vector<std::complex<double> > > &dn);
 
 //helper functions
+template <typename FloatType>
+double eval_delta(const int steps, const double from_value, const double to_value);
 
 
 template<class T> inline T pow2(const T value) {return value*value;}
@@ -70,59 +72,59 @@ int newround(FloatType x) {
   //return x >= 0 ? (x + 0.5).convert_to<int>():(x - 0.5).convert_to<int>();
 }
 template<typename T>
-inline std::complex<T> my_exp(const std::complex<T>& x) {
+inline std::complex<T> my_exp(const std::complex<T> &x) {
   using std::exp; // use ADL
-  T const& r = exp(x.real());
+  T const &r = exp(x.real());
   return std::polar(r, x.imag());
 }
 
 // pl, nmax, mode_n, mode_type
     int nMie(const unsigned int L,
            const int pl,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2,
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2,
            int mode_n, int mode_type);
   // pl and nmax
     int nMie(const unsigned int L,
            const int pl,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
   // no pl and nmax
   int nMie(const unsigned int L,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
   // pl
   int nMie(const unsigned int L,
            const int pl,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
   // nmax
   int nMie(const unsigned int L,
-           std::vector<double>& x, std::vector<std::complex<double> >& m,
-           const unsigned int nTheta, std::vector<double>& Theta,
+           std::vector<double> &x, std::vector<std::complex<double> > &m,
+           const unsigned int nTheta, std::vector<double> &Theta,
            const int nmax,
            double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
            double *g, double *Albedo,
-           std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
+           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
   int nField(const unsigned int L, const int pl,
-             const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax,
+             const std::vector<double> &x, const std::vector<std::complex<double> > &m, const int nmax,
              const int mode_n, const int mode_type, const unsigned int ncoord,
-             const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp,
-             std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
+             const std::vector<double> &Xp, const std::vector<double> &Yp, const std::vector<double> &Zp,
+             std::vector<std::vector<std::complex<double> > > &E, std::vector<std::vector<std::complex<double> > > &H);
 
   // constants for per mode evaluation
   enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};
@@ -142,7 +144,8 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
                                   const int radius_points=1,
                                   const double from_Rho=0, const double to_Rho=static_cast<double>(1.),
                                   const double from_Theta=0, const double to_Theta=static_cast<double>(3.14159265358979323),
-                                  const double from_Phi=0, const double to_Phi=static_cast<double>(3.14159265358979323));
+                                  const double from_Phi=0, const double to_Phi=static_cast<double>(3.14159265358979323),
+                                  const bool isIgnoreAvailableNmax = true); // TODO change to false for production
 
     void calcScattCoeffs();
     void calcExpanCoeffs();
@@ -168,13 +171,15 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
 
     // Problem definition
     // Modify size of all layers
-    void SetLayersSize(const std::vector<FloatType>& layer_size);
+    void SetLayersSize(const std::vector<FloatType> &layer_size);
     // Modify refractive index of all layers
-    void SetLayersIndex(const std::vector< std::complex<FloatType> >& index);
+    void SetLayersIndex(const std::vector< std::complex<FloatType> > &index);
+    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml, unsigned int &l);
+    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml);
     // Modify scattering (theta) angles
-    void SetAngles(const std::vector<FloatType>& angles);
+    void SetAngles(const std::vector<FloatType> &angles);
     // Modify coordinates for field calculation
-    void SetFieldCoords(const std::vector< std::vector<FloatType> >& coords);
+    void SetFieldCoords(const std::vector< std::vector<FloatType> > &coords);
     // Modify index of PEC layer
     void SetPECLayer(int layer_position = 0);
     // Modify the mode taking into account for evaluation of output variables
@@ -236,22 +241,22 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
     std::complex<FloatType> calc_S2(int n, std::complex<FloatType> an, std::complex<FloatType> bn,
                                  FloatType Pi, FloatType Tau);
     void calcD1D3(std::complex<FloatType> z,
-                  std::vector<std::complex<FloatType> >& D1,
-                  std::vector<std::complex<FloatType> >& D3);
+                  std::vector<std::complex<FloatType> > &D1,
+                  std::vector<std::complex<FloatType> > &D3);
     void calcPsiZeta(std::complex<FloatType> x,
-                     std::vector<std::complex<FloatType> >& Psi,
-                     std::vector<std::complex<FloatType> >& Zeta);
-    void calcPiTau(const FloatType& costheta,
-                   std::vector<FloatType>& Pi, std::vector<FloatType>& Tau);
+                     std::vector<std::complex<FloatType> > &Psi,
+                     std::vector<std::complex<FloatType> > &Zeta);
+    void calcPiTau(const FloatType &costheta,
+                   std::vector<FloatType> &Pi, std::vector<FloatType> &Tau);
     void calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
-                       const std::complex<FloatType>& rn, const std::complex<FloatType>& Dn,
-                       const FloatType& Pi, const FloatType& Tau, const FloatType& n,
-                       std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n,
-                       std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n);
+                       const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
+                       const FloatType &Pi, const FloatType &Tau, const FloatType &n,
+                       std::vector<std::complex<FloatType> > &Mo1n, std::vector<std::complex<FloatType> > &Me1n,
+                       std::vector<std::complex<FloatType> > &No1n, std::vector<std::complex<FloatType> > &Ne1n);
 
     void calcFieldByComponents(const FloatType Rho, const FloatType Theta, const FloatType Phi,
-                               std::vector<std::complex<FloatType> >& E,
-                               std::vector<std::complex<FloatType> >& H);
+                               std::vector<std::complex<FloatType> > &E,
+                               std::vector<std::complex<FloatType> > &H);
 
     bool isExpCoeffsCalc_ = false;
     bool isScaCoeffsCalc_ = false;
@@ -268,11 +273,24 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
     // with calcNmax(int first_layer);
     int nmax_ = -1;
     int nmax_preset_ = -1;
+    int available_maximal_nmax_ = -1;
     /// Store result
     FloatType Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
     std::vector<std::vector< std::complex<FloatType> > > E_, H_;  // {X[], Y[], Z[]}
     std::vector<std::vector< std::complex<FloatType> > > Es_, Hs_;  // {X[], Y[], Z[]}
     std::vector<std::complex<FloatType> > S1_, S2_;
+    void calcMieSeriesNeededToConverge(const FloatType Rho);
+    void calcPiTauAllTheta(const double from_Theta,
+                           const double to_Theta,
+                           std::vector<std::vector<FloatType>> &Pi,
+                           std::vector<std::vector<FloatType>> &Tau);
+    void calcRadialOnlyDependantFunctions(const FloatType from_Rho,
+                                          const FloatType to_Rho,
+                                          const bool isIgnoreAvailableNmax,
+                                          std::vector<std::vector<std::complex<FloatType>>> &Psi,
+                                          std::vector<std::vector<std::complex<FloatType>>> &D1n,
+                                          std::vector<std::vector<std::complex<FloatType>>> &Zeta,
+                                          std::vector<std::vector<std::complex<FloatType>>> &D3n);
   };  // end of class MultiLayerMie
 
 }  // end of namespace nmie

src/shell-generator.cc

@@ -45,12 +45,12 @@ namespace shell_generator {
   // ********************************************************************** //
   // ********************************************************************** //
   template<typename T, typename A> inline
-  T dot(std::vector<T,A> const& a,
-        std::vector<T,A> const& b) {
+  T dot(std::vector<T,A> const &a,
+        std::vector<T,A> const &b) {
     return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
     // return std::inner_product(begin(a), end(a), begin(b),
     //                           static_cast<T>(0.0));
-  }  
+  }
   template<typename T, typename A> inline
   std::vector<T> cross(std::vector<T,A> const& a,
         std::vector<T,A> const& b) {
@@ -60,7 +60,7 @@ namespace shell_generator {
       a[0]*b[1]-a[1]*b[0]
     };
     return c;
-  }  
+  }
   // template<typename T, typename A> inline
   // std::vector<std::vector<T> > dyadic(std::vector<T,A> const& a,
   //       std::vector<T,A> const& b) {
@@ -70,7 +70,7 @@ namespace shell_generator {
   //     {a[2]*b[0], a[2]*b[1], a[2]*b[2]}
   //   };
   //   return c;
-  // }  
+  // }
   // template<typename T, typename A> inline
   // std::vector<std::vector<T>,A >
   // operator+(std::vector<std::vector<T>,A > const& a,
@@ -81,7 +81,7 @@ namespace shell_generator {
   //     {a[2][0]+b[2][0], a[2][1]+b[2][1], a[2][2]+b[2][2]}
   //   };
   //   return c;
-  // }  
+  // }
   template<typename T, typename A> inline
   std::vector<T,A >
   operator+(std::vector<T,A > const& a,
@@ -90,7 +90,7 @@ namespace shell_generator {
                            a[1]+b[1],
                            a[2]+b[2]};
     return c;
-  }  
+  }
   template<typename T, typename A> inline
   std::vector<T,A >
   operator-(std::vector<T,A > const& a,
@@ -99,7 +99,7 @@ namespace shell_generator {
                            a[1]-b[1],
                            a[2]-b[2]};
     return c;
-  }  
+  }
 
   // template<typename T, typename A> inline
   // std::vector<std::vector<T>,A >
@@ -111,7 +111,7 @@ namespace shell_generator {
   //     {a[2][0]-b[2][0], a[2][1]-b[2][1], a[2][2]-b[2][2]}
   //   };
   //   return c;
-  // }  
+  // }
   // template<typename T, typename A> inline
   // std::vector<std::vector<T>,A >
   // real(std::vector<std::vector<T>,A > const& a) {
@@ -121,7 +121,7 @@ namespace shell_generator {
   //     {a[2][0].real(), a[2][1].real(), a[2][2].real()}
   //   };
   //   return c;
-  // }  
+  // }
   template<typename T, typename A> inline
   std::vector<T>
   real(std::vector<std::complex<T>,A > const& a) {
@@ -129,12 +129,12 @@ namespace shell_generator {
                         a[1].real(),
                         a[2].real()};
     return c;
-  }  
+  }
   template<typename T, typename A> inline
-  T  real(std::complex<T> const& a) {    
+  T  real(std::complex<T> const& a) {
     return a.real();
-  }  
-  template<typename T, typename A> 
+  }
+  template<typename T, typename A>
   std::vector< std::complex<T>,A > vconj(std::vector< std::complex<T>,A > const& a) {
     std::vector< std::complex<T>,A > b = {std::conj(a[0]),
                                           std::conj(a[1]),
@@ -142,8 +142,8 @@ namespace shell_generator {
     // for (auto elem : a)
     //   b.push_back(std::conj(elem));
     return b;
-  }  
-  template<typename T1, typename T2, typename A> 
+  }
+  template<typename T1, typename T2, typename A>
   std::vector<T2,A> operator*(T1 const& a, std::vector< T2,A > const& b) {
     std::vector<T2,A > c = {a*b[0],
                             a*b[1],
@@ -151,8 +151,8 @@ namespace shell_generator {
     // for (auto elem : b)
     //   c.push_back(a*elem);
     return c;
-  }  
-  template<typename T1, typename T2, typename A> 
+  }
+  template<typename T1, typename T2, typename A>
   std::vector<T2,A> operator/(std::vector< T2,A > const& b, T1 const& a) {
     std::vector<T2,A > c = {b[0]/a,
                             b[1]/a,
@@ -160,10 +160,10 @@ namespace shell_generator {
     // for (auto elem : b)
     //   c.push_back(a*elem);
     return c;
-  }  
+  }
+
+
 
-  
-  
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
@@ -241,7 +241,7 @@ namespace shell_generator {
       std::vector<double> field = {ampl*(p[0]-shift)/r, ampl*(p[1])/r, ampl*(p[2])/r};
       return field;
     };
-    //simple 
+    //simple
     // for (auto vert :vertices_) {
     for (long unsigned int i = 0; i<face_centers_.size(); ++i) {
       auto vert = face_centers_[i];
@@ -253,7 +253,7 @@ namespace shell_generator {
       double r = norm(vert);
       std::vector<double> unit = { vert[0]/r, vert[1]/r, vert[2]/r};
       // std::cout << norm(unit) << std::endl;
-      for (int j =0; j < 3; ++j) 
+      for (int j =0; j < 3; ++j)
         integral += per_face_area_[i]*unit[j]*E0[j];
     }
     return integral;
@@ -269,11 +269,11 @@ namespace shell_generator {
     auto field = [](double charge, double shift, std::vector<double> p){
       double r = std::sqrt(pow2(p[0]-shift) + pow2(p[1]) + pow2(p[2]) );
       const double pi = 3.1415926535897932384626433832795;
-      double ampl = charge/(4.0*pi*pow2(r));      
+      double ampl = charge/(4.0*pi*pow2(r));
       std::vector<double> field = {ampl*(p[0]-shift)/r, ampl*(p[1])/r, ampl*(p[2])/r};
       return field;
     };
-    
+
     for(const auto face : faces_){
       std::vector<double> mean_vector = {0.0, 0.0, 0.0}, mean_point = {0.0, 0.0, 0.0};
       //Get mean
@@ -289,7 +289,7 @@ namespace shell_generator {
       double r = norm(mean_point);
       std::vector<double> unit = { mean_point[0]/r, mean_point[1]/r, mean_point[2]/r};
       // std::cout << norm(unit) << std::endl;
-      for (int i =0; i < 3; ++i) 
+      for (int i =0; i < 3; ++i)
         integral += face_area_*
           unit[i]*mean_vector[i];
     }
@@ -299,7 +299,7 @@ namespace shell_generator {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  
+
   std::vector<double> ShellGenerator::IntegrateByComp() {
     std::vector<double> integral = {0.0, 0.0, 0.0};
     for (unsigned int i=0; i<E_.size(); ++i) {
@@ -311,7 +311,7 @@ namespace shell_generator {
       const std::vector< std::vector<double> > d = {{1.0, 0.0, 0.0},
                                                     {0.0, 1.0, 0.0},
                                                     {0.0, 0.0, 1.0}};
-      
+
       std::vector<double> F = {0.0, 0.0, 0.0};
       std::complex<double> S(0.0);
       for (int ii = 0; ii < 3; ++ii)
@@ -327,14 +327,14 @@ namespace shell_generator {
           F[i] += (1/(2.0/* *4.0*pi */))*real(T[i][j]*n[j]);
         }
       }
-      integral = integral + per_face_area_[i]*F;      
+      integral = integral + per_face_area_[i]*F;
     }
     return integral;
   }
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  
+
   std::vector<double> ShellGenerator::IntegrateByCompReal() {
     std::vector<double> integral = {0.0, 0.0, 0.0};
     for (unsigned int i=0; i<E_.size(); ++i) {
@@ -346,7 +346,7 @@ namespace shell_generator {
       const std::vector< std::vector<double> > d = {{1.0, 0.0, 0.0},
                                                     {0.0, 1.0, 0.0},
                                                     {0.0, 0.0, 1.0}};
-      
+
       std::vector<double> F = {0.0, 0.0, 0.0};
       std::complex<double> S(0.0);
       for (int ii = 0; ii < 3; ++ii)
@@ -362,17 +362,17 @@ namespace shell_generator {
           F[i] += (1/(2.0/* *4.0*pi */))*real(T[i][j]*n[j]);
         }
       }
-      integral = integral + per_face_area_[i]*F;      
+      integral = integral + per_face_area_[i]*F;
     }
     return integral;
   }
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  
+
   std::vector<double> ShellGenerator::IntegrateByFaces() {
     std::vector<double> integral = {0.0, 0.0, 0.0};
-    //simple 
+    //simple
     for (long unsigned int i=0; i<E_.size(); ++i) {
       //std::cout << i << " ";
       auto E = E_[i];
@@ -380,7 +380,7 @@ namespace shell_generator {
       auto H = H_[i];
       // auto Es = Es_[i];
       // auto Hs = Hs_[i];
-      
+
       auto vert = face_centers_[i];
       // Vector to unit product
       double r = norm(vert);
@@ -411,7 +411,7 @@ namespace shell_generator {
       //   (-1.0/4.0)*(dot(E,vconj(E))*unit
       //               +dot(H,vconj(H))*unit
       //               )
-              
+
       //     );
       // auto
       // std::cout <<"E "<<E[0]<<", "<< E[1] <<", "<<E[2] << std::endl;
@@ -434,13 +434,13 @@ namespace shell_generator {
     }
     return integral;
   }
-  
+
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
   std::vector<double> ShellGenerator::Integrate() {
     std::vector<double> integral = {0.0, 0.0, 0.0};
-    //simple 
+    //simple
     for (unsigned int i=0; i<E_.size(); ++i) {
       auto E = E_[i];
       //auto H = 377.0*H_[i];
@@ -471,7 +471,7 @@ namespace shell_generator {
   // ********************************************************************** //
   std::vector< std::vector<double> > ShellGenerator::GetVerticesT() {
     std::vector< std::vector<double> > vertices_t;
-    vertices_t.resize(3); 
+    vertices_t.resize(3);
     for(const auto vert : vertices_){
       vertices_t[0].push_back(vert[0]);
       vertices_t[1].push_back(vert[1]);
@@ -485,7 +485,7 @@ namespace shell_generator {
   std::vector< std::vector<double> > ShellGenerator::GetFaceCentersT() {
     EvalFaces();
     std::vector< std::vector<double> > vertices_t;
-    vertices_t.resize(3); 
+    vertices_t.resize(3);
     for(const auto vert : face_centers_){
       vertices_t[0].push_back(vert[0]);
       vertices_t[1].push_back(vert[1]);
@@ -515,7 +515,7 @@ namespace shell_generator {
        //std::cout << " " << norm(vert) << " ";
      }
      const double pi = 3.1415926535897932384626433832795;
-     double area = 4.0*pi*pow2(scale); 
+     double area = 4.0*pi*pow2(scale);
      //face_area_ = area/faces_.size();
      per_vertice_area_ = area/vertices_.size();
      //std::cout << "Per verice area: " << per_vertice_area_ << std::endl;
@@ -528,7 +528,7 @@ namespace shell_generator {
     for(auto vert : vertices_){
       std::cout <<"(";
       for (auto coord:vert) std::cout<<coord<<",";
-      std::cout <<"),";      
+      std::cout <<"),";
     }
     std::cout << std::endl;
    }
@@ -563,13 +563,13 @@ namespace shell_generator {
       total_flat_area += face;
     auto scale = norm(vertices_[0]);
     const double pi = 3.1415926535897932384626433832795;
-    double area = 4.0*pi*pow2(scale); 
+    double area = 4.0*pi*pow2(scale);
     face_area_ = area/faces_.size();
     double area_scale = area/total_flat_area;
     for (auto& face:per_face_area_)
       face *= area_scale;
-    
-    
+
+
     for(auto& vert : face_centers_){
        double factor = norm(vert);
        //std::cout<< factor <<std::endl;
@@ -634,14 +634,14 @@ namespace shell_generator {
       auto edge_c_index = refined_edges_.size()-1;
 
       /*
-      //                    /\      contrcloсkwise                         
-      //                   c  1                                    
-      //                  0    b                                   
-      //                 /__a___\   edge_a                         
-      //                /\      /\                                 
-      //               c  \    /  0                                
-      //              1    \  /    b                               
-      //    edge_0a  /__0a__\/__1a__\ edge_1a                      
+      //                    /\      contrcloсkwise
+      //                   c  1
+      //                  0    b
+      //                 /__a___\   edge_a
+      //                /\      /\
+      //               c  \    /  0
+      //              1    \  /    b
+      //    edge_0a  /__0a__\/__1a__\ edge_1a
       //
       // remember! In edge_0a the refined point is [1], etc.
       */
@@ -660,7 +660,7 @@ namespace shell_generator {
       // Orient:
       // Try contrcloсkwise:
       bool isClockwise = false, is_b_swapped = false, is_c_swapped=false;
-      
+
       if (edge_0a[0]!=edge_1c[0]) {
         edge_1c.swap(edge_0c);
         is_c_swapped = !is_c_swapped;
@@ -693,18 +693,18 @@ namespace shell_generator {
       }
 
       /*
-      //                    /\      clockwise                               
-      //                   b  1                                    
-      //                  0    c                                   
-      //                 /__a___\   edge_a                         
-      //                /\      /\                                 
-      //               b  \    /  0                                
-      //              1    \  /    c                               
-      //    edge_0a  /__0a__\/__1a__\ edge_1a                      
+      //                    /\      clockwise
+      //                   b  1
+      //                  0    c
+      //                 /__a___\   edge_a
+      //                /\      /\
+      //               b  \    /  0
+      //              1    \  /    c
+      //    edge_0a  /__0a__\/__1a__\ edge_1a
       //
       */
       //Build new facets:
-      // if isClockwise 
+      // if isClockwise
       std::vector<long unsigned int> face1({edge_0a_index, edge_1b_index, edge_c_index});
       std::vector<long unsigned int> face2({edge_1a_index, edge_0c_index, edge_b_index});
       std::vector<long unsigned int> face3({edge_0b_index, edge_1c_index, edge_a_index});
@@ -732,7 +732,7 @@ namespace shell_generator {
       //                << " " << refined_edges_[face1[2]][1] << std::endl;
 
     } // end for faces_
-    
+
     std::cout << "new edges: " << refined_edges_.size() <<std::endl;
     std::cout << "new faces: " << refined_faces_.size() <<std::endl;
     edges_.clear();
@@ -743,7 +743,7 @@ namespace shell_generator {
     // }
     faces_.clear();
     faces_ = refined_faces_;
-    
+
     //Rescale(1.0);
     //GenerateEdges();
     // GenerateFaces();
@@ -766,7 +766,7 @@ namespace shell_generator {
           std::vector<long unsigned int> face({i,j,k});
           // std::cout << ie[0]<<"-"<<ie[1] << ":"
           //             << je[0]<<"-"<<je[1] << ":"
-          //             << ke[0]<<"-"<<ke[1] 
+          //             << ke[0]<<"-"<<ke[1]
           //             << std::endl;
           // std::cout << face[0]<<"-"<<face[1] << "-"<<face[2]<<std::endl;
           faces_.push_back(face);
@@ -837,7 +837,7 @@ namespace shell_generator {
       {a, b,-c},
       {a,-b, c},
       {a,-b,-c},
-      
+
       { b, c,a},
       { b,-c,a},
       {-b, c,a},
@@ -877,7 +877,7 @@ namespace shell_generator {
       {a, b, c},
       {a, -b,-c},
       {-a,-b, c},
-      {-a, b,-c}      
+      {-a, b,-c}
     };
     vertices_ = std::move(points);
     //Rescale(1.0);

src/shell-generator.hpp

@@ -38,9 +38,9 @@
 namespace shell_generator {
   class ShellGenerator {  // will throw for any error
    public:
-    /* ShellGenerator& ReadFromFile(std::string filename); */
-    /* ShellGenerator& ResizeToComplex(double from_wl, double to_wl, int samples); */
-    /* ShellGenerator& ToIndex(); */
+    /* ShellGenerator &ReadFromFile(std::string filename); */
+    /* ShellGenerator &ResizeToComplex(double from_wl, double to_wl, int samples); */
+    /* ShellGenerator &ToIndex(); */
     double dist(std::vector<double> a, std::vector<double> b);
     double norm(std::vector<double> a);
     std::vector< std::vector<double> > GetVertices(){return vertices_;};
@@ -82,7 +82,7 @@ namespace shell_generator {
     std::vector<double> per_face_area_;
     double per_vertice_area_ = 0.0;
     std::vector< std::vector<double> > vertices_, face_centers_;
-    std::vector< double > face_surface_; 
+    std::vector< double > face_surface_;
     std::vector< std::vector<long unsigned int> > edges_, refined_edges_;
     std::vector< std::vector<long unsigned int> > faces_, refined_faces_;
     // std::vector< std::pair< double, std::complex<double> > > data_complex_;

src/special-functions-impl.hpp

@@ -110,7 +110,7 @@ std::complex<FloatType> complex_cot(const std::complex<FloatType> z) {
 // Forward iteration for evaluation of ratio of the Riccati–Bessel functions
 template <typename FloatType>
 void evalForwardR (const std::complex<FloatType> z,
-                   std::vector<std::complex<FloatType> >& r) {
+                   std::vector<std::complex<FloatType> > &r) {
   if (r.size() < 1) throw std::invalid_argument(
       "We need non-zero evaluations of ratio of the Riccati–Bessel functions.\n");
   // r0 = cot(z)
@@ -125,7 +125,7 @@ void evalForwardR (const std::complex<FloatType> z,
 // Backward iteration for evaluation of ratio of the Riccati–Bessel functions
 template <typename FloatType>
 void evalBackwardR (const std::complex<FloatType> z,
-                   std::vector<std::complex<FloatType> >& r) {
+                   std::vector<std::complex<FloatType> > &r) {
   if (r.size() < 1) throw std::invalid_argument(
         "We need non-zero evaluations of ratio of the Riccati–Bessel functions.\n");
   int nmax = r.size()-1 ;
@@ -140,8 +140,8 @@ void evalBackwardR (const std::complex<FloatType> z,
 
 template <typename FloatType>
 void convertRtoD1(const std::complex<FloatType> z,
-                  std::vector<std::complex<FloatType> >& r,
-                  std::vector<std::complex<FloatType> >& D1) {
+                  std::vector<std::complex<FloatType> > &r,
+                  std::vector<std::complex<FloatType> > &D1) {
   if (D1.size() > r.size()) throw std::invalid_argument(
       "Not enough elements in array of ratio of the Riccati–Bessel functions to convert it into logarithmic derivative array.\n");
   std::complex<FloatType> Dold;
@@ -155,7 +155,7 @@ void convertRtoD1(const std::complex<FloatType> z,
 // ********************************************************************** //
 template <typename FloatType>
 void evalForwardD (const std::complex<FloatType> z,
-                     std::vector<std::complex<FloatType> >& D) {
+                     std::vector<std::complex<FloatType> > &D) {
   int nmax = D.size();
   FloatType one = static_cast<FloatType >(1);
   for (int n = 1; n < nmax; n++) {
@@ -167,7 +167,7 @@ void evalForwardD (const std::complex<FloatType> z,
 // ********************************************************************** //
 template <typename FloatType>
 void evalForwardD1 (const std::complex<FloatType> z,
-                   std::vector<std::complex<FloatType> >& D) {
+                   std::vector<std::complex<FloatType> > &D) {
   if (D.size()<1) throw std::invalid_argument("Should have a leas one element!\n");
   D[0] = std::cos(z)/std::sin(z);
   evalForwardD(z,D);
@@ -175,8 +175,8 @@ void evalForwardD1 (const std::complex<FloatType> z,
 
 //  template <typename FloatType>
 //  void MultiLayerMie<FloatType>::calcD1D3(const std::complex<FloatType> z,
-//                               std::vector<std::complex<FloatType> >& D1,
-//                               std::vector<std::complex<FloatType> >& D3) {
+//                               std::vector<std::complex<FloatType> > &D1,
+//                               std::vector<std::complex<FloatType> > &D3) {
 //
 //    // if (cabs(D1[0]) > 1.0e15) {
 //    //   throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
@@ -210,8 +210,8 @@ void evalForwardD1 (const std::complex<FloatType> z,
   //**********************************************************************************//
 //  template <typename FloatType>
 //  void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
-//                                  std::vector<std::complex<FloatType> >& Psi,
-//                                  std::vector<std::complex<FloatType> >& Zeta) {
+//                                  std::vector<std::complex<FloatType> > &Psi,
+//                                  std::vector<std::complex<FloatType> > &Zeta) {
 //
 //    std::complex<FloatType> c_i(0.0, 1.0);
 //    std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1);
@@ -243,8 +243,8 @@ void evalForwardD1 (const std::complex<FloatType> z,
   //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
   //**********************************************************************************//
 //  template <typename FloatType>
-//  void MultiLayerMie<FloatType>::calcPiTau(const FloatType& costheta,
-//                                std::vector<FloatType>& Pi, std::vector<FloatType>& Tau) {
+//  void MultiLayerMie<FloatType>::calcPiTau(const FloatType &costheta,
+//                                std::vector<FloatType> &Pi, std::vector<FloatType> &Tau) {
 //
 //    int i;
 //    //****************************************************//
@@ -283,8 +283,8 @@ void evalForwardD1 (const std::complex<FloatType> z,
   //**********************************************************************************//
 //  template <typename FloatType>
 //  void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
-//                                    const std::complex<FloatType>& rn, const std::complex<FloatType>& Dn,
-//                                    const FloatType& Pi, const FloatType& Tau, const FloatType& n,
+//                                    const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
+//                                    const FloatType &Pi, const FloatType &Tau, const FloatType &n,
 //                                    std::vector<std::complex<FloatType> >& Mo1n, std::vector<std::complex<FloatType> >& Me1n,
 //                                    std::vector<std::complex<FloatType> >& No1n, std::vector<std::complex<FloatType> >& Ne1n) {
 //

tests/c++/speed-test-applied.cc

@@ -83,19 +83,19 @@ int main(int argc, char *argv[]) {
     double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow;
 
     double ti = 0.0, tf = 90.0;
-    int nt = 0;    
+    int nt = 0;
     if (argc < 5) throw std::invalid_argument(error_msg);
-    
+
     //strcpy(comment, "");
     // for (i = 1; i < argc; i++) {
-    int mode = -1; 
+    int mode = -1;
     double tmp_mr;
     for (auto arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
 
-      // Detecting new read mode (if it is a valid -key) 
+      // Detecting new read mode (if it is a valid -key)
       if (arg == "-l") {
 	mode = read_L;
 	continue;
@@ -160,13 +160,13 @@ int main(int argc, char *argv[]) {
 	continue;
       }
     }
-    if ( (x.size() != m.size()) || (L != x.size()) ) 
+    if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n")
 							 +error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) ) 
+    if ( (0 == m.size()) || ( 0 == x.size()) )
       throw std::invalid_argument(std::string("Empty structure!\n")
 							 +error_msg);
-    
+
     if (nt < 0) {
       printf("Error reading Theta.\n");
       return -1;
@@ -182,7 +182,7 @@ int main(int argc, char *argv[]) {
     long ctime_nsec, best_c;
     long cpptime_sec, ctime_sec;
     long repeats = 150;
-    //HeapProfilerStart("heapprof");    
+    //HeapProfilerStart("heapprof");
     do {
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
       for (int i = 0; i<repeats; ++i) {
@@ -198,26 +198,26 @@ int main(int argc, char *argv[]) {
     } while (cpptime_nsec < 1e8 && ctime_nsec < 1e8);
 
         printf("\n");
-    
+
     if (has_comment) {
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     } else {
       printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     }
-    
+
     if (nt > 0) {
       printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
-      
+
       for (i = 0; i < nt; i++) {
         printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
       }
     }
 
-  } 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;
     return -1;
-  }  
+  }
     return 0;
 }
 

tests/c++/speed-test.cc

@@ -85,19 +85,19 @@ int main(int argc, char *argv[]) {
     double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow;
 
     double ti = 0.0, tf = 90.0;
-    int nt = 0;    
+    int nt = 0;
     if (argc < 5) throw std::invalid_argument(error_msg);
-    
+
     //strcpy(comment, "");
     // for (i = 1; i < argc; i++) {
-    int mode = -1; 
+    int mode = -1;
     double tmp_mr;
     for (auto arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
 
-      // Detecting new read mode (if it is a valid -key) 
+      // Detecting new read mode (if it is a valid -key)
       if (arg == "-l") {
 	mode = read_L;
 	continue;
@@ -162,13 +162,13 @@ int main(int argc, char *argv[]) {
 	continue;
       }
     }
-    if ( (x.size() != m.size()) || (L != x.size()) ) 
+    if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n")
 							 +error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) ) 
+    if ( (0 == m.size()) || ( 0 == x.size()) )
       throw std::invalid_argument(std::string("Empty structure!\n")
 							 +error_msg);
-    
+
     if (nt < 0) {
       printf("Error reading Theta.\n");
       return -1;
@@ -184,7 +184,7 @@ int main(int argc, char *argv[]) {
     long ctime_nsec, best_c;
     long cpptime_sec, ctime_sec;
     long repeats = 150;
-    //HeapProfilerStart("heapprof");    
+    //HeapProfilerStart("heapprof");
     printf("Start\n");
     do {
       clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
@@ -203,26 +203,26 @@ int main(int argc, char *argv[]) {
     } while (cpptime_sec < 3 && ctime_sec < 3);
 
         printf("Finish\n");
-    
+
     // if (has_comment) {
     //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     // } else {
     //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     // }
-    
+
     // if (nt > 0) {
     //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
-      
+
     //   for (i = 0; i < nt; i++) {
     //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
     //   }
     // }
 
-  } 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;
     return -1;
-  }  
+  }
     return 0;
 }
 

tests/test_near_field.cc

@@ -18,7 +18,7 @@ TEST(RunFieldCalculationPolar, HandlesInput) {
   nmie.SetLayersSize({0.099});
   nmie.SetLayersIndex({ {0.75,0}});
   nmie.RunMieCalculation();
-  nmie.RunFieldCalculationPolar(2);
+//  nmie.RunFieldCalculationPolar(2);
 }
 //TEST(BulkSphere, HandlesInput) {
 //  nmie::MultiLayerMie<nmie::FloatType> nmie;

utils/bessel/bessel.cc

@@ -35,8 +35,8 @@
 namespace nmie {
   namespace bessel {
 
-    void calcZeta(int n,  std::complex<double>z,  std::vector< std::complex<double> >& Zeta,
-		   std::vector< std::complex<double> >& dZeta) {
+    void calcZeta(int n,  std::complex<double>z,  std::vector< std::complex<double> > &Zeta,
+		   std::vector< std::complex<double> > &dZeta) {
       std::vector< std::complex<double> > csj, cdj, csy, cdy;
       int nm;
       csphjy (n, z, nm, csj, cdj,  csy, cdy );
@@ -49,8 +49,8 @@ namespace nmie {
       }
     }  // end of calcZeta()
 
-    void calcPsi(int n,  std::complex<double>z,  std::vector< std::complex<double> >& Psi,
-		   std::vector< std::complex<double> >& dPsi) {
+    void calcPsi(int n,  std::complex<double>z,  std::vector< std::complex<double> > &Psi,
+		   std::vector< std::complex<double> > &dPsi) {
       std::vector< std::complex<double> > csj, cdj, csy, cdy;
       int nm;
       csphjy (n, z, nm, csj, cdj,  csy, cdy );
@@ -63,9 +63,9 @@ namespace nmie {
     }  // end of calcPsi()
 
 // !*****************************************************************************80
-//  
+//
 //  C++ port of fortran code
-//  
+//
 // !! CSPHJY: spherical Bessel functions jn(z) and yn(z) for complex argument.
 // !
 // !  Discussion:
@@ -75,8 +75,8 @@ namespace nmie {
 // !
 // !  Licensing:
 // !
-// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
-// !    they give permission to incorporate this routine into a user program 
+// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+// !    they give permission to incorporate this routine into a user program
 // !    provided that the copyright is acknowledged.
 // !
 // !  Modified:
@@ -94,7 +94,7 @@ namespace nmie {
 // !    Wiley, 1996,
 // !    ISBN: 0-471-11963-6,
 // !    LC: QA351.C45.
-// ! 
+// !
 // !  Parameters:
 // !
 // !    Input, integer ( kind = 4 ) N, the order of jn(z) and yn(z).
@@ -106,11 +106,11 @@ namespace nmie {
 // !    Output, complex ( kind = 8 ) CSJ(0:N0, CDJ(0:N), CSY(0:N), CDY(0:N),
 // !    the values of jn(z), jn'(z), yn(z), yn'(z).
 // !
-    void csphjy (int n, std::complex<double>z, int& nm,
-		 std::vector< std::complex<double> >& csj,
-		 std::vector< std::complex<double> >& cdj,
-		 std::vector< std::complex<double> >& csy,
-		 std::vector< std::complex<double> >& cdy ) {
+    void csphjy (int n, std::complex<double>z, int &nm,
+		 std::vector< std::complex<double> > &csj,
+		 std::vector< std::complex<double> > &cdj,
+		 std::vector< std::complex<double> > &csy,
+		 std::vector< std::complex<double> > &cdy ) {
       double a0;
       csj.resize(n+1);
       cdj.resize(n+1);
@@ -119,7 +119,7 @@ namespace nmie {
       std::complex<double> cf, cf0, cf1, cs, csa, csb;
       int m;
       a0 = std::abs(z);
-      nm = n;      
+      nm = n;
       if (a0 < 1.0e-60) {
 	for (int k = 0; k < n+1; ++k) {
 	  csj[k] = 0.0;
@@ -133,7 +133,7 @@ namespace nmie {
       }
       csj[0] = std::sin ( z ) / z;
       csj[1] = ( csj[0] - std::cos ( z ) ) / z;
-      
+
       if ( 2 <= n ) {
 	csa = csj[0];
 	csb = csj[1];
@@ -174,11 +174,11 @@ namespace nmie {
       for (int k = 2; k<=nm; ++k) {
 	cdy[k] = csy[k-1] - ( k + 1.0 ) * csy[k] / z;
       }
-      
+
       return;
     }
       // function msta2 ( x, n, mp )
-      
+
       // !*****************************************************************************80
       // !
       // !! MSTA2 determines a backward recurrence starting point for Jn(x).
@@ -190,8 +190,8 @@ namespace nmie {
       // !
       // !  Licensing:
       // !
-      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
-      // !    they give permission to incorporate this routine into a user program 
+      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+      // !    they give permission to incorporate this routine into a user program
       // !    provided that the copyright is acknowledged.
       // !
       // !  Modified:
@@ -257,14 +257,14 @@ namespace nmie {
 // !
 // !  Discussion:
 // !
-// !    This procedure determines the starting point for backward  
-// !    recurrence such that the magnitude of    
+// !    This procedure determines the starting point for backward
+// !    recurrence such that the magnitude of
 // !    Jn(x) at that point is about 10^(-MP).
 // !
 // !  Licensing:
 // !
-// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
-// !    they give permission to incorporate this routine into a user program 
+// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+// !    they give permission to incorporate this routine into a user program
 // !    provided that the copyright is acknowledged.
 // !
 // !  Modified:
@@ -287,7 +287,7 @@ namespace nmie {
 // !
 // !    Input, real ( kind = 8 ) X, the argument.
 // !
-// !    Input, integer ( kind = 4 ) MP, the negative logarithm of the 
+// !    Input, integer ( kind = 4 ) MP, the negative logarithm of the
 // !    desired magnitude.
 // !
 // !    Output, integer ( kind = 4 ) MSTA1, the starting point.
@@ -300,7 +300,7 @@ namespace nmie {
       f0 = envj ( n0, a0 ) - mp;
       n1 = n0 + 5;
       f1 = envj ( n1, a0 ) - mp;
-      for (int it = 1; it <= 20; ++it) { 
+      for (int it = 1; it <= 20; ++it) {
 	nn = n1 - ( n1 - n0 ) / ( 1.0 - f0 / f1 );
 	f = envj ( nn, a0 ) - mp;
 	if ( abs ( nn - n1 ) < 1 ) break;
@@ -319,8 +319,8 @@ namespace nmie {
       // !
       // !  Licensing:
       // !
-      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
-      // !    they give permission to incorporate this routine into a user program 
+      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+      // !    they give permission to incorporate this routine into a user program
       // !    provided that the copyright is acknowledged.
       // !
       // !  Modified:

utils/bessel/bessel.h

@@ -29,8 +29,8 @@
 namespace nmie {
   namespace bessel {
     void calcZeta(int n,  std::complex<double>z,
-		  std::vector< std::complex<double> >& Zeta,
-		  std::vector< std::complex<double> >& dZeta);
+		  std::vector< std::complex<double> > &Zeta,
+		  std::vector< std::complex<double> > &dZeta);
     void calcPsi(int n,  std::complex<double>z,
 		  std::vector< std::complex<double> >& Psi,
 		  std::vector< std::complex<double> >& dPsi);

utils/bessel/test_bessel.cc

@@ -37,7 +37,7 @@ int nmax_ = -1;
 //Temporary variables
 std::vector<std::complex<double> > PsiZeta_;
 
-void calcD1D3(const std::complex<double> z, std::vector<std::complex<double> >& D1, std::vector<std::complex<double> >& D3) {
+void calcD1D3(const std::complex<double> z, std::vector<std::complex<double> > &D1, std::vector<std::complex<double> > &D3) {
     D1[nmax_] = std::complex<double>(0.0, 0.0);
     const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
     for (int n = nmax_; n > 0; n--) {
@@ -59,7 +59,7 @@ void calcD1D3(const std::complex<double> z, std::vector<std::complex<double> >&
     }
   }
 
-  void calcPsiZeta(std::complex<double> z, std::vector<std::complex<double> >& Psi, std::vector<std::complex<double> >& Zeta) {
+  void calcPsiZeta(std::complex<double> z, std::vector<std::complex<double> > &Psi, std::vector<std::complex<double> > &Zeta) {
     std::complex<double> c_i(0.0, 1.0);
     std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
     calcD1D3(z, D1, D3);
@@ -76,7 +76,7 @@ void calcD1D3(const std::complex<double> z, std::vector<std::complex<double> >&
 
 
 int main(int argc, char *argv[]) {
-  
+
   int n_small = 1, n_big = 15;
   int n = n_big+2, n_print;
   std::complex<double> z_small(0.0012, 0.0034);
@@ -104,7 +104,7 @@ int main(int argc, char *argv[]) {
   for (unsigned int i = 0; i < csj.size(); ++i)
     csh[i] = csj[i]+c_i*csy[i];
   nmie::bessel::calcPsi(n, z, Psi, dPsi);
-  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);
 
   // n_print = n_small;
   // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
@@ -125,7 +125,7 @@ int main(int argc, char *argv[]) {
   // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 csh[n_print].real(), csh[n_print].imag());
   // printf("WA csj  =         e   r        e   i;\n");
-    
+
   n_print = n_small;
   printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
   	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
@@ -174,7 +174,7 @@ int main(int argc, char *argv[]) {
   for (unsigned int i = 0; i < csj.size(); ++i)
     csh[i] = csj[i]+c_i*csy[i];
   nmie::bessel::calcPsi(n, z, Psi, dPsi);
-  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);
 
   // n_print = n_small;
   // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
@@ -195,7 +195,7 @@ int main(int argc, char *argv[]) {
   // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 csh[n_print].real(), csh[n_print].imag());
   // printf("WA csj  =         e   r        e   i;\n");
-    
+
   n_print = n_small;
   printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
   	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
@@ -241,7 +241,7 @@ int main(int argc, char *argv[]) {
   for (unsigned int i = 0; i < csj.size(); ++i)
     csh[i] = csj[i]+c_i*csy[i];
   nmie::bessel::calcPsi(n, z, Psi, dPsi);
-  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);
 
   // n_print = n_small;
   // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
@@ -262,7 +262,7 @@ int main(int argc, char *argv[]) {
   // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 csh[n_print].real(), csh[n_print].imag());
   // printf("WA csj  =         e   r        e   i;\n");
-    
+
   n_print = n_small;
   printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
   	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
@@ -309,7 +309,7 @@ int main(int argc, char *argv[]) {
   for (unsigned int i = 0; i < csj.size(); ++i)
     csh[i] = csj[i]+c_i*csy[i];
   nmie::bessel::calcPsi(n, z, Psi, dPsi);
-  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);
 
   // n_print = n_small;
   // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
@@ -330,7 +330,7 @@ int main(int argc, char *argv[]) {
   // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 csh[n_print].real(), csh[n_print].imag());
   // printf("WA csj  =         e   r        e   i;\n");
-    
+
   n_print = n_small;
   printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
   	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
@@ -380,7 +380,7 @@ int main(int argc, char *argv[]) {
   for (unsigned int i = 0; i < csj.size(); ++i)
     csh[i] = csj[i]+c_i*csy[i];
   nmie::bessel::calcPsi(n, z, Psi, dPsi);
-  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);
 
   // n_print = n_small;
   // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
@@ -401,7 +401,7 @@ int main(int argc, char *argv[]) {
   // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 csh[n_print].real(), csh[n_print].imag());
   // printf("WA csj  =         e   r        e   i;\n");
-    
+
   // n_print = n_small;
   // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
   // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,

utils/compare.cc

@@ -45,9 +45,9 @@ const double PI=3.14159265358979323846;
 template<class T> inline T pow2(const T value) {return value*value;}
 
 template <class VectorType, int dimensions> inline
-std::vector<VectorType> CrossProduct(std::vector<VectorType>& a, std::vector<VectorType>& b) {
+std::vector<VectorType> CrossProduct(std::vector<VectorType> &a, std::vector<VectorType> &b) {
   if (a.size() != 3 || b.size() != 3) throw std::invalid_argument("Cross product only for 3D vectors!");
-  std::vector<VectorType> r (3);   
+  std::vector<VectorType> r (3);
   r[0] = a[1]*b[2]-a[2]*b[1];
   r[1] = a[2]*b[0]-a[0]*b[2];
   r[2] = a[0]*b[1]-a[1]*b[0];
@@ -92,19 +92,19 @@ int main(int argc, char *argv[]) {
     double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow;
 
     double ti = 0.0, tf = 90.0;
-    int nt = 0;    
+    int nt = 0;
     if (argc < 5) throw std::invalid_argument(error_msg);
-    
+
     //strcpy(comment, "");
     // for (i = 1; i < argc; i++) {
-    int mode = -1; 
+    int mode = -1;
     double tmp_mr;
     for (auto arg : args) {
       // For each arg in args list we detect the change of the current
       // read mode or read the arg. The reading args algorithm works
       // as a finite-state machine.
 
-      // Detecting new read mode (if it is a valid -key) 
+      // Detecting new read mode (if it is a valid -key)
       if (arg == "-l") {
 	mode = read_L;
 	continue;
@@ -169,13 +169,13 @@ int main(int argc, char *argv[]) {
 	continue;
       }
     }
-    if ( (x.size() != m.size()) || (L != x.size()) ) 
+    if ( (x.size() != m.size()) || (L != x.size()) )
       throw std::invalid_argument(std::string("Broken structure!\n")
 							 +error_msg);
-    if ( (0 == m.size()) || ( 0 == x.size()) ) 
+    if ( (0 == m.size()) || ( 0 == x.size()) )
       throw std::invalid_argument(std::string("Empty structure!\n")
 							 +error_msg);
-    
+
     if (nt < 0) {
       printf("Error reading Theta.\n");
       return -1;
@@ -191,7 +191,7 @@ int main(int argc, char *argv[]) {
     // long ctime_nsec, best_c;
     // long cpptime_sec, ctime_sec;
     // long repeats = 150;
-    // //HeapProfilerStart("heapprof");    
+    // //HeapProfilerStart("heapprof");
     // do {
     //   clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
     //   for (int i = 0; i<repeats; ++i) {
@@ -202,9 +202,9 @@ int main(int argc, char *argv[]) {
     //   cpptime_nsec = diff(time1,time2).tv_nsec;
     //   cpptime_sec = diff(time1,time2).tv_sec;
     //   clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
-    //   // for (int i = 0; i<repeats; ++i) {      
+    //   // for (int i = 0; i<repeats; ++i) {
     //   // 	nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-    //   // }  
+    //   // }
     //   clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
     //   ctime_nsec = diff(time1,time2).tv_nsec;
     //   ctime_sec = diff(time1,time2).tv_sec;
@@ -226,7 +226,7 @@ int main(int argc, char *argv[]) {
     nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
     nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
         printf("\n");
-    
+
     if (has_comment) {
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  old\n", comment.c_str(), Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
@@ -234,10 +234,10 @@ int main(int argc, char *argv[]) {
       printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  old\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
       printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     }
-    
+
     if (nt > 0) {
       printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
-      
+
       for (i = 0; i < nt; i++) {
         printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  old\n", Theta[i]*180.0/PI, S1[i].real(), S1[i].imag(), S2[i].real(), S2[i].imag());
         printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
@@ -297,8 +297,8 @@ int main(int argc, char *argv[]) {
     int pl = 0;
     int nmax = 0;
     std::vector<std::vector<std::complex<double> > > E(ncoord), H(ncoord);
-    for (auto& f:E) f.resize(3);
-    for (auto& f:H) f.resize(3);
+    for (auto &f:E) f.resize(3);
+    for (auto &f:H) f.resize(3);
     double free_impedance = 376.73031;
     //double free_impedance = 1.0;
 
@@ -320,7 +320,7 @@ int main(int argc, char *argv[]) {
       }
       if (sum_e > max_E) max_E = sum_e;
       if (sum_e < min_E) min_E = sum_e;
-    
+
       //printf("Field E=%g\n", std::sqrt(std::abs(sum_e)));
     }
     printf("Min E = %g; max E =%g", min_E, max_E);
@@ -330,11 +330,11 @@ int main(int argc, char *argv[]) {
     //   printf("Field H=%g\n", std::sqrt(std::abs(sum_h))*free_impedance);
     // }
 
-  } 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;
     return -1;
-  }  
+  }
     return 0;
 }
 

utils/test-negative-epsilon.cc

@@ -2,9 +2,9 @@
  * @file   test-negative-epsilon.cc
  * @author Konstantin Ladutenko <kostyfisik at gmail (.) com>
  * @date   Mon Mar  9 13:21:37 2015
- * 
+ *
  * @brief  test negative epsilon case
- * 
+ *
 //    This file is part of scattnlay                                                //
 //                                                                                  //
 //    This program is free software: you can redistribute it and/or modify          //
@@ -26,13 +26,13 @@
 //                                                                                  //
 //    You should have received a copy of the GNU General Public License             //
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
- * 
+ *
  */
 #include "nmie.h"
 #include <stdio.h>
 template<class T> inline T pow2(const T value) {return value*value;};
-const double PI=3.14159265358979323846;  
-nmie::MultiLayerMie multi_layer_mie_;  
+const double PI=3.14159265358979323846;
+nmie::MultiLayerMie multi_layer_mie_;
 double lambda_work_ = 3.75; // cm
 double a_ = 0.75*lambda_work_;  // 2.8125 cm - size of PEC core
 double min_index_ = 1e-11;
@@ -59,12 +59,12 @@ double EvaluateScatterOnlyIndex(std::vector<double> input) {
   try {
     multi_layer_mie_.RunMieCalculations();
     Qsca = multi_layer_mie_.GetQsca();
-  } catch( const std::invalid_argument& ia ) {
+  } catch( const std::invalid_argument &ia ) {
     Qsca = 0;
     printf("#");
     // Will catch if  multi_layer_mie_ fails or other errors.
     //std::cerr << "Invalid argument: " << ia.what() << std::endl;
-  }  
+  }
   double total_r = multi_layer_mie_.GetTotalRadius();
   return Qsca*PI*pow2(total_r);
 }
@@ -78,7 +78,7 @@ int main(int argc, char *argv[]) {
     double PEC_Qsca = multi_layer_mie_.GetQsca();
     double PEC_r = multi_layer_mie_.GetTotalRadius();
     double PEC_RCS = PEC_Qsca*PI*pow2(PEC_r);
-    
+
     // PEC target covered with with air layer
     multi_layer_mie_.SetCoatingWidth({0.1});
     multi_layer_mie_.SetCoatingIndex({{1.0,0.0}});
@@ -87,7 +87,7 @@ int main(int argc, char *argv[]) {
     double total_r1 = multi_layer_mie_.GetTotalRadius();
     double initial_RCS1 = Qsca1*PI*pow2(total_r1);
     printf("RCS = %g cm^2 with (r=%g) and  RCS=%g cm^2 without (r=%g)air coating.\n",
-	   initial_RCS1, total_r1, 
+	   initial_RCS1, total_r1,
 	   PEC_RCS, PEC_r);
 
     //multi_layer_mie.SetMaxTermsNumber(150);
@@ -131,13 +131,13 @@ int main(int argc, char *argv[]) {
 	   EvaluateScatterOnlyIndex({-0.29, 24.6, 1.0}));
     //multi_layer_mie_.SetMaxTermsNumber(-1);
 
-      // 26.24: 25||   -0.29   +24.60 
+      // 26.24: 25||   -0.29   +24.60
       // 28.48: 38||   -0.29   +24.60    +1.00
-    
-  } 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;
     return -1;
-  }  
+  }
     return 0;
 }