|  | @@ -45,6 +45,8 @@
 | 
	
		
			
				|  |  |  #include <stdexcept>
 | 
	
		
			
				|  |  |  #include <vector>
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | +using namespace std;
 | 
	
		
			
				|  |  | +
 | 
	
		
			
				|  |  |  namespace nmie {
 | 
	
		
			
				|  |  |    //helpers
 | 
	
		
			
				|  |  |    template<class T> inline T pow2(const T value) {return value*value;}
 | 
	
	
		
			
				|  | @@ -72,10 +74,10 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Return value:                                                                    //
 | 
	
		
			
				|  |  |    //   Number of multipolar expansion terms used for the calculations                 //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  | -  int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
 | 
	
		
			
				|  |  | +  int ScattCoeffs(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const int nmax, vector<complex<double> >& an, vector<complex<double> >& bn) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      if (x.size() != L || m.size() != L)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Declared number of layers do not fit x and m!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Declared number of layers do not fit x and m!");
 | 
	
		
			
				|  |  |      try {
 | 
	
		
			
				|  |  |        MultiLayerMie ml_mie;
 | 
	
		
			
				|  |  |        ml_mie.SetAllLayersSize(x);
 | 
	
	
		
			
				|  | @@ -89,10 +91,10 @@ namespace nmie {
 | 
	
		
			
				|  |  |        bn = ml_mie.GetBn();
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        return ml_mie.GetMaxTerms();
 | 
	
		
			
				|  |  | -    } catch(const std::invalid_argument& ia) {
 | 
	
		
			
				|  |  | +    } catch(const 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);
 | 
	
		
			
				|  |  | +      cerr << "Invalid argument: " << ia.what() << endl;
 | 
	
		
			
				|  |  | +      throw invalid_argument(ia);
 | 
	
		
			
				|  |  |        return -1;
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |      return 0;
 | 
	
	
		
			
				|  | @@ -127,12 +129,12 @@ 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, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
 | 
	
		
			
				|  |  | +  int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      if (x.size() != L || m.size() != L)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Declared number of layers do not fit x and m!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Declared number of layers do not fit x and m!");
 | 
	
		
			
				|  |  |      if (Theta.size() != nTheta)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Declared number of sample for Theta is not correct!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Declared number of sample for Theta is not correct!");
 | 
	
		
			
				|  |  |      try {
 | 
	
		
			
				|  |  |        MultiLayerMie ml_mie;
 | 
	
		
			
				|  |  |        ml_mie.SetAllLayersSize(x);
 | 
	
	
		
			
				|  | @@ -154,10 +156,10 @@ namespace nmie {
 | 
	
		
			
				|  |  |        S2 = ml_mie.GetS2();
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        return ml_mie.GetMaxTerms();
 | 
	
		
			
				|  |  | -    } catch(const std::invalid_argument& ia) {
 | 
	
		
			
				|  |  | +    } catch(const 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);
 | 
	
		
			
				|  |  | +      cerr << "Invalid argument: " << ia.what() << endl;
 | 
	
		
			
				|  |  | +      throw invalid_argument(ia);
 | 
	
		
			
				|  |  |        return -1;
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |      return 0;
 | 
	
	
		
			
				|  | @@ -191,7 +193,7 @@ 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, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
 | 
	
		
			
				|  |  | +  int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
 | 
	
		
			
				|  |  |      return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -223,7 +225,7 @@ 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, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
 | 
	
		
			
				|  |  | +  int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
 | 
	
		
			
				|  |  |      return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -257,7 +259,7 @@ 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, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
 | 
	
		
			
				|  |  | +  int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
 | 
	
		
			
				|  |  |      return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -284,18 +286,18 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Return value:                                                                    //
 | 
	
		
			
				|  |  |    //   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 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) {
 | 
	
		
			
				|  |  | +  int nField(const unsigned int L, const int pl, const vector<double>& x, const vector<complex<double> >& m, const int nmax, const unsigned int ncoord, const vector<double>& Xp_vec, const vector<double>& Yp_vec, const vector<double>& Zp_vec, vector<vector<complex<double> > >& E, vector<vector<complex<double> > >& H) {
 | 
	
		
			
				|  |  |      if (x.size() != L || m.size() != L)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Declared number of layers do not fit x and m!");
 | 
	
		
			
				|  |  | +      throw 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!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
 | 
	
		
			
				|  |  |      for (auto f:E)
 | 
	
		
			
				|  |  |        if (f.size() != 3)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Field E is not 3D!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Field E is not 3D!");
 | 
	
		
			
				|  |  |      for (auto f:H)
 | 
	
		
			
				|  |  |        if (f.size() != 3)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Field H is not 3D!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Field H is not 3D!");
 | 
	
		
			
				|  |  |      try {
 | 
	
		
			
				|  |  |        MultiLayerMie ml_mie;
 | 
	
		
			
				|  |  |        //ml_mie.SetPECLayer(pl); // TODO add PEC layer to field plotting
 | 
	
	
		
			
				|  | @@ -307,10 +309,10 @@ namespace nmie {
 | 
	
		
			
				|  |  |        H = ml_mie.GetFieldH();
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        return ml_mie.GetMaxTerms();
 | 
	
		
			
				|  |  | -    } catch(const std::invalid_argument& ia) {
 | 
	
		
			
				|  |  | +    } catch(const 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);
 | 
	
		
			
				|  |  | +      cerr << "Invalid argument: " << ia.what() << endl;
 | 
	
		
			
				|  |  | +      throw invalid_argument(ia);
 | 
	
		
			
				|  |  |        return - 1;
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |      return 0;
 | 
	
	
		
			
				|  | @@ -322,7 +324,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetQext() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return Qext_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -332,7 +334,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetQabs() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return Qabs_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -342,7 +344,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetQsca() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return Qsca_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -352,7 +354,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetQbk() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return Qbk_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -362,7 +364,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetQpr() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return Qpr_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -372,7 +374,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetAsymmetryFactor() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return asymmetry_factor_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -382,7 +384,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    double MultiLayerMie::GetAlbedo() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return albedo_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -390,9 +392,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Returns previously calculated S1                                       //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::vector<std::complex<double> > MultiLayerMie::GetS1() {
 | 
	
		
			
				|  |  | +  vector<complex<double> > MultiLayerMie::GetS1() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return S1_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -400,9 +402,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Returns previously calculated S2                                       //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::vector<std::complex<double> > MultiLayerMie::GetS2() {
 | 
	
		
			
				|  |  | +  vector<complex<double> > MultiLayerMie::GetS2() {
 | 
	
		
			
				|  |  |      if (!isMieCalculated_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("You should run calculations before result request!");
 | 
	
		
			
				|  |  |      return S2_;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -410,7 +412,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Modify scattering (theta) angles                                       //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  void MultiLayerMie::SetAngles(const std::vector<double>& angles) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::SetAngles(const vector<double>& angles) {
 | 
	
		
			
				|  |  |      isExpCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isMieCalculated_ = false;
 | 
	
	
		
			
				|  | @@ -421,7 +423,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Modify size of all layers                                             //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  void MultiLayerMie::SetAllLayersSize(const std::vector<double>& layer_size) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::SetAllLayersSize(const vector<double>& layer_size) {
 | 
	
		
			
				|  |  |      isExpCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isMieCalculated_ = false;
 | 
	
	
		
			
				|  | @@ -429,9 +431,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |      double prev_layer_size = 0.0;
 | 
	
		
			
				|  |  |      for (auto curr_layer_size : layer_size) {
 | 
	
		
			
				|  |  |        if (curr_layer_size <= 0.0)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument("Size parameter should be positive!");
 | 
	
		
			
				|  |  | +        throw invalid_argument("Size parameter should be positive!");
 | 
	
		
			
				|  |  |        if (prev_layer_size > curr_layer_size)
 | 
	
		
			
				|  |  | -        throw std::invalid_argument
 | 
	
		
			
				|  |  | +        throw invalid_argument
 | 
	
		
			
				|  |  |            ("Size parameter for next layer should be larger than the previous one!");
 | 
	
		
			
				|  |  |        prev_layer_size = curr_layer_size;
 | 
	
		
			
				|  |  |        size_param_.push_back(curr_layer_size);
 | 
	
	
		
			
				|  | @@ -442,7 +444,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Modify refractive index of all layers                                  //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  void MultiLayerMie::SetAllLayersIndex(const std::vector< std::complex<double> >& index) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::SetAllLayersIndex(const vector< complex<double> >& index) {
 | 
	
		
			
				|  |  |      isExpCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isMieCalculated_ = false;
 | 
	
	
		
			
				|  | @@ -453,11 +455,11 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Modify coordinates for field calculation                               //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  void MultiLayerMie::SetFieldCoords(const std::vector< std::vector<double> >& coords) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::SetFieldCoords(const vector< vector<double> >& coords) {
 | 
	
		
			
				|  |  |      if (coords.size() != 3)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Error! Wrong dimension of field monitor points!");
 | 
	
		
			
				|  |  |      if (coords[0].size() != coords[1].size() || coords[0].size() != coords[2].size())
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Error! Missing coordinates for field monitor points!");
 | 
	
		
			
				|  |  |      coords_ = coords;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -470,7 +472,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isMieCalculated_ = false;
 | 
	
		
			
				|  |  |      if (layer_position < 0 && layer_position != -1)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Error! Layers are numbered from 0!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Error! Layers are numbered from 0!");
 | 
	
		
			
				|  |  |      PEC_layer_position_ = layer_position;
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -538,21 +540,21 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    void MultiLayerMie::calcNmax(unsigned int first_layer) {
 | 
	
		
			
				|  |  |      int ri, riM1;
 | 
	
		
			
				|  |  | -    const std::vector<double>& x = size_param_;
 | 
	
		
			
				|  |  | -    const std::vector<std::complex<double> >& m = refractive_index_;
 | 
	
		
			
				|  |  | +    const vector<double>& x = size_param_;
 | 
	
		
			
				|  |  | +    const vector<complex<double> >& m = refractive_index_;
 | 
	
		
			
				|  |  |      calcNstop();  // Set initial nmax_ value
 | 
	
		
			
				|  |  |      for (unsigned int i = first_layer; i < x.size(); i++) {
 | 
	
		
			
				|  |  |        if (static_cast<int>(i) > PEC_layer_position_)  // static_cast used to avoid warning
 | 
	
		
			
				|  |  | -        ri = round(std::abs(x[i]*m[i]));
 | 
	
		
			
				|  |  | +        ri = round(abs(x[i]*m[i]));
 | 
	
		
			
				|  |  |        else
 | 
	
		
			
				|  |  |          ri = 0;
 | 
	
		
			
				|  |  | -      nmax_ = std::max(nmax_, ri);
 | 
	
		
			
				|  |  | +      nmax_ = max(nmax_, ri);
 | 
	
		
			
				|  |  |        // first layer is pec, if pec is present
 | 
	
		
			
				|  |  |        if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
 | 
	
		
			
				|  |  | -        riM1 = round(std::abs(x[i - 1]* m[i]));
 | 
	
		
			
				|  |  | +        riM1 = round(abs(x[i - 1]* m[i]));
 | 
	
		
			
				|  |  |        else
 | 
	
		
			
				|  |  |          riM1 = 0;
 | 
	
		
			
				|  |  | -      nmax_ = std::max(nmax_, riM1);
 | 
	
		
			
				|  |  | +      nmax_ = max(nmax_, riM1);
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |      nmax_ += 15;  // Final nmax_ value
 | 
	
		
			
				|  |  |    }
 | 
	
	
		
			
				|  | @@ -561,12 +563,12 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Calculate an - equation (5)                                            //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::complex<double> MultiLayerMie::calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
 | 
	
		
			
				|  |  | -                                              std::complex<double> PsiXL, std::complex<double> ZetaXL,
 | 
	
		
			
				|  |  | -                                              std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
 | 
	
		
			
				|  |  | +  complex<double> MultiLayerMie::calc_an(int n, double XL, complex<double> Ha, complex<double> mL,
 | 
	
		
			
				|  |  | +                                         complex<double> PsiXL, complex<double> ZetaXL,
 | 
	
		
			
				|  |  | +                                         complex<double> PsiXLM1, complex<double> ZetaXLM1) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
 | 
	
		
			
				|  |  | -    std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
 | 
	
		
			
				|  |  | +    complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
 | 
	
		
			
				|  |  | +    complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      return Num/Denom;
 | 
	
		
			
				|  |  |    }
 | 
	
	
		
			
				|  | @@ -575,12 +577,12 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Calculate bn - equation (6)                                            //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::complex<double> MultiLayerMie::calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
 | 
	
		
			
				|  |  | -                                              std::complex<double> PsiXL, std::complex<double> ZetaXL,
 | 
	
		
			
				|  |  | -                                              std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
 | 
	
		
			
				|  |  | +  complex<double> MultiLayerMie::calc_bn(int n, double XL, complex<double> Hb, complex<double> mL,
 | 
	
		
			
				|  |  | +                                         complex<double> PsiXL, complex<double> ZetaXL,
 | 
	
		
			
				|  |  | +                                         complex<double> PsiXLM1, complex<double> ZetaXLM1) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
 | 
	
		
			
				|  |  | -    std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
 | 
	
		
			
				|  |  | +    complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
 | 
	
		
			
				|  |  | +    complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      return Num/Denom;
 | 
	
		
			
				|  |  |    }
 | 
	
	
		
			
				|  | @@ -589,8 +591,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  |    // Calculates S1 - equation (25a)                                         //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::complex<double> MultiLayerMie::calc_S1(int n, std::complex<double> an, std::complex<double> bn,
 | 
	
		
			
				|  |  | -                                              double Pi, double Tau) {
 | 
	
		
			
				|  |  | +  complex<double> MultiLayerMie::calc_S1(int n, complex<double> an, complex<double> bn,
 | 
	
		
			
				|  |  | +                                         double Pi, double Tau) {
 | 
	
		
			
				|  |  | +
 | 
	
		
			
				|  |  |      return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -599,8 +602,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Calculates S2 - equation (25b) (it's the same as (25a), just switches  //
 | 
	
		
			
				|  |  |    // Pi and Tau)                                                            //
 | 
	
		
			
				|  |  |    // ********************************************************************** //
 | 
	
		
			
				|  |  | -  std::complex<double> MultiLayerMie::calc_S2(int n, std::complex<double> an, std::complex<double> bn,
 | 
	
		
			
				|  |  | -                                              double Pi, double Tau) {
 | 
	
		
			
				|  |  | +  complex<double> MultiLayerMie::calc_S2(int n, complex<double> an, complex<double> bn,
 | 
	
		
			
				|  |  | +                                         double Pi, double Tau) {
 | 
	
		
			
				|  |  | +
 | 
	
		
			
				|  |  |      return calc_S1(n, an, bn, Tau, Pi);
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -617,31 +621,31 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Output parameters:                                                               //
 | 
	
		
			
				|  |  |    //   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions                //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  | -  void MultiLayerMie::calcD1D3(const std::complex<double> z,
 | 
	
		
			
				|  |  | -                               std::vector<std::complex<double> >& D1,
 | 
	
		
			
				|  |  | -                               std::vector<std::complex<double> >& D3) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::calcD1D3(const complex<double> z,
 | 
	
		
			
				|  |  | +                               vector<complex<double> >& D1,
 | 
	
		
			
				|  |  | +                               vector<complex<double> >& D3) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Downward recurrence for D1 - equations (16a) and (16b)
 | 
	
		
			
				|  |  | -    D1[nmax_] = std::complex<double>(0.0, 0.0);
 | 
	
		
			
				|  |  | -    const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
 | 
	
		
			
				|  |  | +    D1[nmax_] = complex<double>(0.0, 0.0);
 | 
	
		
			
				|  |  | +    const complex<double> zinv = complex<double>(1.0, 0.0)/z;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      for (int n = nmax_; n > 0; n--) {
 | 
	
		
			
				|  |  |        D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    if (std::abs(D1[0]) > 1.0e15) {
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
 | 
	
		
			
				|  |  | +    if (abs(D1[0]) > 1.0e15) {
 | 
	
		
			
				|  |  | +      throw invalid_argument("Unstable D1! Please, try to change input parameters!\n");
 | 
	
		
			
				|  |  |      //printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
 | 
	
		
			
				|  |  | -    PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
 | 
	
		
			
				|  |  | -                      *std::exp(-2.0*z.imag()));
 | 
	
		
			
				|  |  | -    D3[0] = std::complex<double>(0.0, 1.0);
 | 
	
		
			
				|  |  | +    PsiZeta_[0] = 0.5*(1.0 - complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))
 | 
	
		
			
				|  |  | +                      *exp(-2.0*z.imag()));
 | 
	
		
			
				|  |  | +    D3[0] = complex<double>(0.0, 1.0);
 | 
	
		
			
				|  |  |      for (int n = 1; n <= nmax_; n++) {
 | 
	
		
			
				|  |  |        PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
 | 
	
		
			
				|  |  |                                     *(static_cast<double>(n)*zinv - D3[n - 1]);
 | 
	
		
			
				|  |  | -      D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
 | 
	
		
			
				|  |  | +      D3[n] = D1[n] + complex<double>(0.0, 1.0)/PsiZeta_[n];
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |    }
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -658,19 +662,19 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Output parameters:                                                               //
 | 
	
		
			
				|  |  |    //   Psi, Zeta: Riccati-Bessel functions                                            //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  | -  void MultiLayerMie::calcPsiZeta(std::complex<double> z,
 | 
	
		
			
				|  |  | -                                  std::vector<std::complex<double> >& Psi,
 | 
	
		
			
				|  |  | -                                  std::vector<std::complex<double> >& Zeta) {
 | 
	
		
			
				|  |  | +  void MultiLayerMie::calcPsiZeta(complex<double> z,
 | 
	
		
			
				|  |  | +                                  vector<complex<double> >& Psi,
 | 
	
		
			
				|  |  | +                                  vector<complex<double> >& Zeta) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> c_i(0.0, 1.0);
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
 | 
	
		
			
				|  |  | +    complex<double> c_i(0.0, 1.0);
 | 
	
		
			
				|  |  | +    vector<complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // First, calculate the logarithmic derivatives
 | 
	
		
			
				|  |  |      calcD1D3(z, D1, D3);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
 | 
	
		
			
				|  |  | -    Psi[0] = std::sin(z);
 | 
	
		
			
				|  |  | -    Zeta[0] = std::sin(z) - c_i*std::cos(z);
 | 
	
		
			
				|  |  | +    Psi[0] = sin(z);
 | 
	
		
			
				|  |  | +    Zeta[0] = sin(z) - c_i*cos(z);
 | 
	
		
			
				|  |  |      for (int n = 1; n <= nmax_; n++) {
 | 
	
		
			
				|  |  |        Psi[n]  =  Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
 | 
	
		
			
				|  |  |        Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
 | 
	
	
		
			
				|  | @@ -692,7 +696,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |    //   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c)   //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  |    void MultiLayerMie::calcPiTau(const double& costheta,
 | 
	
		
			
				|  |  | -                                std::vector<double>& Pi, std::vector<double>& Tau) {
 | 
	
		
			
				|  |  | +                                vector<double>& Pi, vector<double>& Tau) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      int i;
 | 
	
		
			
				|  |  |      //****************************************************//
 | 
	
	
		
			
				|  | @@ -729,17 +733,15 @@ namespace nmie {
 | 
	
		
			
				|  |  |    // Output parameters:                                                               //
 | 
	
		
			
				|  |  |    //   Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics                     //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  | -  void MultiLayerMie::calcSpherHarm(const std::complex<double> Rho, const double Theta, const double Phi,
 | 
	
		
			
				|  |  | -                                    const std::complex<double>& rn, const std::complex<double>& Dn,
 | 
	
		
			
				|  |  | +  void MultiLayerMie::calcSpherHarm(const complex<double> Rho, const double Theta, const double Phi,
 | 
	
		
			
				|  |  | +                                    const complex<double>& rn, const complex<double>& Dn,
 | 
	
		
			
				|  |  |                                      const double& Pi, const double& Tau, const double& n,
 | 
	
		
			
				|  |  | -                                    std::vector<std::complex<double> >& Mo1n, std::vector<std::complex<double> >& Me1n, 
 | 
	
		
			
				|  |  | -                                    std::vector<std::complex<double> >& No1n, std::vector<std::complex<double> >& Ne1n) {
 | 
	
		
			
				|  |  | +                                    vector<complex<double> >& Mo1n, vector<complex<double> >& Me1n, 
 | 
	
		
			
				|  |  | +                                    vector<complex<double> >& No1n, vector<complex<double> >& Ne1n) {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // using eq 4.50 in BH
 | 
	
		
			
				|  |  | -    std::complex<double> c_zero(0.0, 0.0);
 | 
	
		
			
				|  |  | +    complex<double> c_zero(0.0, 0.0);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    using std::sin;
 | 
	
		
			
				|  |  | -    using std::cos;
 | 
	
		
			
				|  |  |      Mo1n[0] = c_zero;
 | 
	
		
			
				|  |  |      Mo1n[1] = cos(Phi)*Pi*rn/Rho;
 | 
	
		
			
				|  |  |      Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
 | 
	
	
		
			
				|  | @@ -778,8 +780,8 @@ namespace nmie {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    const std::vector<double>& x = size_param_;
 | 
	
		
			
				|  |  | -    const std::vector<std::complex<double> >& m = refractive_index_;
 | 
	
		
			
				|  |  | +    const vector<double>& x = size_param_;
 | 
	
		
			
				|  |  | +    const vector<complex<double> >& m = refractive_index_;
 | 
	
		
			
				|  |  |      const int& pl = PEC_layer_position_;
 | 
	
		
			
				|  |  |      const int L = refractive_index_.size();
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -792,7 +794,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |      if (nmax_preset_ <= 0) calcNmax(fl);
 | 
	
		
			
				|  |  |      else nmax_ = nmax_preset_;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> z1, z2;
 | 
	
		
			
				|  |  | +    complex<double> z1, z2;
 | 
	
		
			
				|  |  |      //**************************************************************************//
 | 
	
		
			
				|  |  |      // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which  //
 | 
	
		
			
				|  |  |      // means that index = layer number - 1 or index = n - 1. The only exception //
 | 
	
	
		
			
				|  | @@ -801,10 +803,10 @@ namespace nmie {
 | 
	
		
			
				|  |  |      // between different arrays. The change was done to optimize memory usage.  //
 | 
	
		
			
				|  |  |      //**************************************************************************//
 | 
	
		
			
				|  |  |      // Allocate memory to the arrays
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
 | 
	
		
			
				|  |  | +    vector<complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
 | 
	
		
			
				|  |  |                                         D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::vector<std::vector<std::complex<double> > > Q(L), Ha(L), Hb(L);
 | 
	
		
			
				|  |  | +    vector<vector<complex<double> > > Q(L), Ha(L), Hb(L);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      for (int l = 0; l < L; l++) {
 | 
	
		
			
				|  |  |        Q[l].resize(nmax_ + 1);
 | 
	
	
		
			
				|  | @@ -816,15 +818,15 @@ namespace nmie {
 | 
	
		
			
				|  |  |      bn_.resize(nmax_);
 | 
	
		
			
				|  |  |      PsiZeta_.resize(nmax_ + 1);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
 | 
	
		
			
				|  |  | +    vector<complex<double> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      //*************************************************//
 | 
	
		
			
				|  |  |      // Calculate D1 and D3 for z1 in the first layer   //
 | 
	
		
			
				|  |  |      //*************************************************//
 | 
	
		
			
				|  |  |      if (fl == pl) {  // PEC layer
 | 
	
		
			
				|  |  |        for (int n = 0; n <= nmax_; n++) {
 | 
	
		
			
				|  |  | -        D1_mlxl[n] = std::complex<double>(0.0, - 1.0);
 | 
	
		
			
				|  |  | -        D3_mlxl[n] = std::complex<double>(0.0, 1.0);
 | 
	
		
			
				|  |  | +        D1_mlxl[n] = complex<double>(0.0, - 1.0);
 | 
	
		
			
				|  |  | +        D3_mlxl[n] = complex<double>(0.0, 1.0);
 | 
	
		
			
				|  |  |        }
 | 
	
		
			
				|  |  |      } else { // Regular layer
 | 
	
		
			
				|  |  |        z1 = x[fl]* m[fl];
 | 
	
	
		
			
				|  | @@ -842,8 +844,8 @@ namespace nmie {
 | 
	
		
			
				|  |  |      //*****************************************************//
 | 
	
		
			
				|  |  |      // Iteration from the second layer to the last one (L) //
 | 
	
		
			
				|  |  |      //*****************************************************//
 | 
	
		
			
				|  |  | -    std::complex<double> Temp, Num, Denom;
 | 
	
		
			
				|  |  | -    std::complex<double> G1, G2;
 | 
	
		
			
				|  |  | +    complex<double> Temp, Num, Denom;
 | 
	
		
			
				|  |  | +    complex<double> G1, G2;
 | 
	
		
			
				|  |  |      for (int l = fl + 1; l < L; l++) {
 | 
	
		
			
				|  |  |        //************************************************************//
 | 
	
		
			
				|  |  |        //Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L   //
 | 
	
	
		
			
				|  | @@ -859,9 +861,9 @@ namespace nmie {
 | 
	
		
			
				|  |  |        //Calculate Q, Ha and Hb in the layers fl + 1..L   //
 | 
	
		
			
				|  |  |        //*************************************************//
 | 
	
		
			
				|  |  |        // Upward recurrence for Q - equations (19a) and (19b)
 | 
	
		
			
				|  |  | -      Num = std::exp(-2.0*(z1.imag() - z2.imag()))
 | 
	
		
			
				|  |  | -           *std::complex<double>(std::cos(-2.0*z2.real()) - std::exp(-2.0*z2.imag()), std::sin(-2.0*z2.real()));
 | 
	
		
			
				|  |  | -      Denom = std::complex<double>(std::cos(-2.0*z1.real()) - std::exp(-2.0*z1.imag()), std::sin(-2.0*z1.real()));
 | 
	
		
			
				|  |  | +      Num = exp(-2.0*(z1.imag() - z2.imag()))
 | 
	
		
			
				|  |  | +           *complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
 | 
	
		
			
				|  |  | +      Denom = complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
 | 
	
		
			
				|  |  |        Q[l][0] = Num/Denom;
 | 
	
		
			
				|  |  |        for (int n = 1; n <= nmax_; n++) {
 | 
	
		
			
				|  |  |          Num = (z1*D1_mlxl[n] + double(n))*(double(n) - z1*D3_mlxl[n - 1]);
 | 
	
	
		
			
				|  | @@ -919,7 +921,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |          an_[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
 | 
	
		
			
				|  |  |          bn_[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
 | 
	
		
			
				|  |  |        } else {
 | 
	
		
			
				|  |  | -        an_[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
 | 
	
		
			
				|  |  | +        an_[n] = calc_an(n + 1, x[L - 1], complex<double>(0.0, 0.0), complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
 | 
	
		
			
				|  |  |          bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
 | 
	
		
			
				|  |  |        }
 | 
	
		
			
				|  |  |      }  // end of for an and bn terms
 | 
	
	
		
			
				|  | @@ -957,11 +959,11 @@ namespace nmie {
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  |    void MultiLayerMie::RunMieCalculation() {
 | 
	
		
			
				|  |  |      if (size_param_.size() != refractive_index_.size())
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Each size parameter should have only one index!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Each size parameter should have only one index!");
 | 
	
		
			
				|  |  |      if (size_param_.size() == 0)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Initialize model first!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Initialize model first!");
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    const std::vector<double>& x = size_param_;
 | 
	
		
			
				|  |  | +    const vector<double>& x = size_param_;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      isExpCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |      isScaCoeffsCalc_ = false;
 | 
	
	
		
			
				|  | @@ -971,7 +973,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |      calcScattCoeffs();
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      if (!isScaCoeffsCalc_) // TODO seems to be unreachable
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("Calculation of scattering coefficients failed!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("Calculation of scattering coefficients failed!");
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Initialize the scattering parameters
 | 
	
		
			
				|  |  |      Qext_ = 0.0;
 | 
	
	
		
			
				|  | @@ -983,14 +985,14 @@ namespace nmie {
 | 
	
		
			
				|  |  |      albedo_ = 0.0;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Initialize the scattering amplitudes
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > tmp1(theta_.size(),std::complex<double>(0.0, 0.0));
 | 
	
		
			
				|  |  | +    vector<complex<double> > tmp1(theta_.size(),complex<double>(0.0, 0.0));
 | 
	
		
			
				|  |  |      S1_.swap(tmp1);
 | 
	
		
			
				|  |  |      S2_ = S1_;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::vector<double> Pi(nmax_), Tau(nmax_);
 | 
	
		
			
				|  |  | +    vector<double> Pi(nmax_), Tau(nmax_);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> Qbktmp(0.0, 0.0);
 | 
	
		
			
				|  |  | -    std::vector< std::complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
 | 
	
		
			
				|  |  | +    complex<double> Qbktmp(0.0, 0.0);
 | 
	
		
			
				|  |  | +    vector< complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
 | 
	
		
			
				|  |  |      // By using downward recurrence we avoid loss of precision due to float rounding errors
 | 
	
		
			
				|  |  |      // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
 | 
	
		
			
				|  |  |      //      http://en.wikipedia.org/wiki/Loss_of_significance
 | 
	
	
		
			
				|  | @@ -1002,14 +1004,14 @@ namespace nmie {
 | 
	
		
			
				|  |  |        Qsca_ += (n + n + 1.0)*(an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
 | 
	
		
			
				|  |  |                              + bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
 | 
	
		
			
				|  |  |        // Equation (29)
 | 
	
		
			
				|  |  | -      Qpr_ += ((n*(n + 2)/(n + 1))*((an_[i]*std::conj(an_[n]) + bn_[i]*std::conj(bn_[n])).real())
 | 
	
		
			
				|  |  | -               + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*std::conj(bn_[i])).real());
 | 
	
		
			
				|  |  | +      Qpr_ += ((n*(n + 2)/(n + 1))*((an_[i]*conj(an_[n]) + bn_[i]*conj(bn_[n])).real())
 | 
	
		
			
				|  |  | +               + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*conj(bn_[i])).real());
 | 
	
		
			
				|  |  |        // Equation (33)
 | 
	
		
			
				|  |  |        Qbktmp += (double)(n + n + 1)*(1 - 2*(n % 2))*(an_[i]- bn_[i]);
 | 
	
		
			
				|  |  |        // Calculate the scattering amplitudes (S1 and S2)    //
 | 
	
		
			
				|  |  |        // Equations (25a) - (25b)                            //
 | 
	
		
			
				|  |  |        for (unsigned int t = 0; t < theta_.size(); t++) {
 | 
	
		
			
				|  |  | -        calcPiTau(std::cos(theta_[t]), Pi, Tau);
 | 
	
		
			
				|  |  | +        calcPiTau(cos(theta_[t]), Pi, Tau);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |          S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
 | 
	
		
			
				|  |  |          S2_[t] += calc_S2(n, an_[i], bn_[i], Pi[i], Tau[i]);
 | 
	
	
		
			
				|  | @@ -1049,11 +1051,11 @@ namespace nmie {
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  |    void MultiLayerMie::calcExpanCoeffs() {
 | 
	
		
			
				|  |  |      if (!isScaCoeffsCalc_)
 | 
	
		
			
				|  |  | -      throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
 | 
	
		
			
				|  |  | +      throw invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      isExpCoeffsCalc_ = false;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
 | 
	
		
			
				|  |  | +    complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      const int L = refractive_index_.size();
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -1077,17 +1079,17 @@ namespace nmie {
 | 
	
		
			
				|  |  |        dln_[L][n] = c_one;
 | 
	
		
			
				|  |  |      }
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
 | 
	
		
			
				|  |  | -    std::complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
 | 
	
		
			
				|  |  | +    vector<complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
 | 
	
		
			
				|  |  | +    vector<complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
 | 
	
		
			
				|  |  | +    complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      auto& m = refractive_index_;
 | 
	
		
			
				|  |  | -    std::vector< std::complex<double> > m1(L);
 | 
	
		
			
				|  |  | +    vector< complex<double> > m1(L);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      for (int l = 0; l < L - 1; l++) m1[l] = m[l + 1];
 | 
	
		
			
				|  |  | -    m1[L - 1] = std::complex<double> (1.0, 0.0);
 | 
	
		
			
				|  |  | +    m1[L - 1] = complex<double> (1.0, 0.0);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> z, z1;
 | 
	
		
			
				|  |  | +    complex<double> z, z1;
 | 
	
		
			
				|  |  |      for (int l = L - 1; l >= 0; l--) {
 | 
	
		
			
				|  |  |        z = size_param_[l]*m[l];
 | 
	
		
			
				|  |  |        z1 = size_param_[l]*m1[l];
 | 
	
	
		
			
				|  | @@ -1122,15 +1124,15 @@ namespace nmie {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Check the result and change  aln_[0][n] and aln_[0][n] for exact zero
 | 
	
		
			
				|  |  |      for (int n = 0; n < nmax_; ++n) {
 | 
	
		
			
				|  |  | -      if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
 | 
	
		
			
				|  |  | +      if (abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
 | 
	
		
			
				|  |  |        else {
 | 
	
		
			
				|  |  | -        //throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
 | 
	
		
			
				|  |  | +        //throw invalid_argument("Unstable calculation of aln_[0][n]!");
 | 
	
		
			
				|  |  |          printf("Warning: Potentially unstable calculation of aln (aln[0][%i] = %g, %gi)\n", n, aln_[0][n].real(), aln_[0][n].imag());
 | 
	
		
			
				|  |  |          aln_[0][n] = 0.0;
 | 
	
		
			
				|  |  |        }
 | 
	
		
			
				|  |  | -      if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
 | 
	
		
			
				|  |  | +      if (abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
 | 
	
		
			
				|  |  |        else {
 | 
	
		
			
				|  |  | -        //throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
 | 
	
		
			
				|  |  | +        //throw invalid_argument("Unstable calculation of bln_[0][n]!");
 | 
	
		
			
				|  |  |          printf("Warning: Potentially unstable calculation of bln (bln[0][%i] = %g, %gi)\n", n, bln_[0][n].real(), bln_[0][n].imag());
 | 
	
		
			
				|  |  |          bln_[0][n] = 0.0;
 | 
	
		
			
				|  |  |        }
 | 
	
	
		
			
				|  | @@ -1153,17 +1155,17 @@ namespace nmie {
 | 
	
		
			
				|  |  |    //   E, H: Complex electric and magnetic fields                                     //
 | 
	
		
			
				|  |  |    //**********************************************************************************//
 | 
	
		
			
				|  |  |    void MultiLayerMie::calcField(const double Rho, const double Theta, const double Phi,
 | 
	
		
			
				|  |  | -                                std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H)  {
 | 
	
		
			
				|  |  | +                                vector<complex<double> >& E, vector<complex<double> >& H)  {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    std::complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
 | 
	
		
			
				|  |  | -    std::vector<std::complex<double> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
 | 
	
		
			
				|  |  | -    std::vector<double> Pi(nmax_), Tau(nmax_);
 | 
	
		
			
				|  |  | +    complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
 | 
	
		
			
				|  |  | +    vector<complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
 | 
	
		
			
				|  |  | +    vector<complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
 | 
	
		
			
				|  |  | +    vector<complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
 | 
	
		
			
				|  |  | +    vector<complex<double> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
 | 
	
		
			
				|  |  | +    vector<double> Pi(nmax_), Tau(nmax_);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      int l = 0;  // Layer number
 | 
	
		
			
				|  |  | -    std::complex<double> ml;
 | 
	
		
			
				|  |  | +    complex<double> ml;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Initialize E and H
 | 
	
		
			
				|  |  |      for (int i = 0; i < 3; i++) {
 | 
	
	
		
			
				|  | @@ -1189,7 +1191,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |      calcPsiZeta(Rho*ml, Psi, Zeta);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // Calculate angular functions Pi and Tau
 | 
	
		
			
				|  |  | -    calcPiTau(std::cos(Theta), Pi, Tau);
 | 
	
		
			
				|  |  | +    calcPiTau(cos(Theta), Pi, Tau);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      for (int n = nmax_ - 2; n >= 0; n--) {
 | 
	
		
			
				|  |  |        int n1 = n + 1;
 | 
	
	
		
			
				|  | @@ -1200,7 +1202,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |        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<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
 | 
	
		
			
				|  |  | +      complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
 | 
	
		
			
				|  |  |        for (int i = 0; i < 3; i++) {
 | 
	
		
			
				|  |  |          // electric field E [V m - 1] = EF*E0
 | 
	
		
			
				|  |  |          E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
 | 
	
	
		
			
				|  | @@ -1212,7 +1214,7 @@ namespace nmie {
 | 
	
		
			
				|  |  |      }  // end of for all n
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      // magnetic field
 | 
	
		
			
				|  |  | -    std::complex<double> hffact = ml/(cc_*mu_);
 | 
	
		
			
				|  |  | +    complex<double> hffact = ml/(cc_*mu_);
 | 
	
		
			
				|  |  |      for (int i = 0; i < 3; i++) {
 | 
	
		
			
				|  |  |        H[i] = hffact*H[i];
 | 
	
		
			
				|  |  |      }
 | 
	
	
		
			
				|  | @@ -1262,16 +1264,16 @@ namespace nmie {
 | 
	
		
			
				|  |  |        const double& Zp = coords_[2][point];
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // Convert to spherical coordinates
 | 
	
		
			
				|  |  | -      Rho = std::sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
 | 
	
		
			
				|  |  | +      Rho = sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
 | 
	
		
			
				|  |  | -      Theta = (Rho > 0.0) ? std::acos(Zp/Rho) : 0.0;
 | 
	
		
			
				|  |  | +      Theta = (Rho > 0.0) ? acos(Zp/Rho) : 0.0;
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // If Xp=Yp=0 then Phi is undefined. Just set it to zero to avoid problems
 | 
	
		
			
				|  |  |        if (Xp == 0.0)
 | 
	
		
			
				|  |  | -        Phi = (Yp != 0.0) ? std::asin(Yp/std::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
 | 
	
		
			
				|  |  | +        Phi = (Yp != 0.0) ? asin(Yp/sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
 | 
	
		
			
				|  |  |        else
 | 
	
		
			
				|  |  | -        Phi = std::acos(Xp/std::sqrt(pow2(Xp) + pow2(Yp)));
 | 
	
		
			
				|  |  | +        Phi = acos(Xp/sqrt(pow2(Xp) + pow2(Yp)));
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // Avoid convergence problems due to Rho too small
 | 
	
		
			
				|  |  |        if (Rho < 1e-5) Rho = 1e-5;
 | 
	
	
		
			
				|  | @@ -1283,14 +1285,12 @@ namespace nmie {
 | 
	
		
			
				|  |  |        //*******************************************************//
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // This array contains the fields in spherical coordinates
 | 
	
		
			
				|  |  | -      std::vector<std::complex<double> > Es(3), Hs(3);
 | 
	
		
			
				|  |  | +      vector<complex<double> > Es(3), Hs(3);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        // Do the actual calculation of electric and magnetic field
 | 
	
		
			
				|  |  |        calcField(Rho, Theta, Phi, Es, Hs);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        { //Now, convert the fields back to cartesian coordinates
 | 
	
		
			
				|  |  | -        using std::sin;
 | 
	
		
			
				|  |  | -        using std::cos;
 | 
	
		
			
				|  |  |          E_[point][0] = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
 | 
	
		
			
				|  |  |          E_[point][1] = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
 | 
	
		
			
				|  |  |          E_[point][2] = cos(Theta)*Es[0] - sin(Theta)*Es[1];
 |