add argument to calcNmax()

src/nmie-basic.hpp

@@ -275,7 +275,7 @@ namespace nmie {
   // ********************************************************************** //
 
   template <typename FloatType>
-  int LeRu_cutoff(const std::complex<FloatType> zz) {
+  unsigned int LeRu_cutoff(const std::complex<FloatType> zz) {
     std::complex<double> z = ConvertComplex<double>(zz);
     auto x = std::abs(z);
     return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
@@ -286,19 +286,21 @@ namespace nmie {
   // Calculate calcNstop - equation (17)                                    //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcNstop() {
+  unsigned int MultiLayerMie<FloatType>::calcNstop(FloatType xL) {
+    unsigned int nmax = 0;
     //Wiscombe
-    const FloatType& xL = size_param_.back();
+    if (xL < size_param_.back()) xL = size_param_.back();
     if (xL <= 8) {
-      nmax_ = newround(xL + 4.0*pow(xL, 1.0/3.0) + 1);
+      nmax = newround(xL + 4.0*pow(xL, 1.0/3.0) + 1);
     } else if (xL <= 4200) {
-      nmax_ = newround(xL + 4.05*pow(xL, 1.0/3.0) + 2);
+      nmax = newround(xL + 4.05*pow(xL, 1.0/3.0) + 2);
     } else {
-      nmax_ = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
+      nmax = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
     }
     //Le Ru
     auto Nstop = nmie::LeRu_cutoff(std::complex<FloatType>(xL,0))+1;
-    if (Nstop > nmax_) nmax_ = Nstop;
+    if (Nstop > nmax) nmax = Nstop;
+    return nmax;
   }
 
 
@@ -306,30 +308,33 @@ namespace nmie {
   // Maximum number of terms required for the calculation                   //
   // ********************************************************************** //
   template <typename FloatType>
-  void MultiLayerMie<FloatType>::calcNmax(unsigned int first_layer) {
-    int ri, riM1;
+  unsigned int MultiLayerMie<FloatType>::calcNmax(FloatType xL) {
+    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_;
-    calcNstop();  // Set initial nmax_ value
+    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
         ri = newround(cabs(x[i]*m[i]));
       else
         ri = 0;
-      nmax_ = std::max(nmax_, ri);
+      nmax = std::max(nmax, ri);
       // first layer is pec, if pec is present
       if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
         riM1 = newround(cabs(x[i - 1]* m[i]));
       else
         riM1 = 0;
-      nmax_ = std::max(nmax_, riM1);
+      nmax = std::max(nmax, riM1);
     }
-    nmax_ += 100;  // Final nmax_ value
+    nmax += 15;  // Final nmax value
 #ifdef MULTI_PRECISION
-    nmax_ += MULTI_PRECISION; //TODO we may need to use more terms that this for MP computations.
+    nmax += MULTI_PRECISION; //TODO we may need to use more terms that this for MP computations.
 #endif
-    // nmax_ *= nmax_;
-    // printf("using nmax %i\n", nmax_);
+    // nmax *= nmax;
+    // printf("using nmax %i\n", nmax);
+    return nmax;
   }
 
 
@@ -459,7 +464,10 @@ namespace nmie {
   void MultiLayerMie<FloatType>::calcPiTau(const FloatType& costheta,
                                 std::vector<FloatType>& Pi, std::vector<FloatType>& Tau) {
 
-    int i;
+    int nmax = Pi.size();
+    if (Pi.size() != Tau.size())
+      throw std::invalid_argument("Error! Pi and Tau vectors should have the same size!");
+
     //****************************************************//
     // Equations (26a) - (26c)                            //
     //****************************************************//
@@ -467,10 +475,10 @@ namespace nmie {
     Pi[0] = 1.0;  // n=1
     Tau[0] = costheta;
     // Calculate the actual values
-    if (nmax_ > 1) {
+    if (nmax > 1) {
       Pi[1] = 3*costheta*Pi[0]; //n=2
       Tau[1] = 2*costheta*Pi[1] - 3*Pi[0];
-      for (i = 2; i < nmax_; i++) { //n=[3..nmax_]
+      for (int i = 2; i < nmax; i++) { //n=[3..nmax_]
         Pi[i] = ((i + i + 1)*costheta*Pi[i - 1] - (i + 1)*Pi[i - 2])/i;
         Tau[i] = (i + 1)*costheta*Pi[i] - (i + 2)*Pi[i - 1];
       }
@@ -560,7 +568,7 @@ namespace nmie {
     // below the PEC are discarded.                                           //
     // ***********************************************************************//
     int fl = (pl > 0) ? pl : 0;
-    if (nmax_preset_ <= 0) calcNmax(fl);
+    if (nmax_preset_ <= 0) nmax_ = calcNmax();
     else nmax_ = nmax_preset_;
 
     std::complex<FloatType> z1, z2;

src/nmie-nearfield.hpp

@@ -425,24 +425,24 @@ double eval_delta(const int steps, const double from_value, const double to_valu
 
 
 // input parameters:
-//         in_outer_perimeter_points: will be increased to the nearest power of 2.
+//         input_outer_perimeter_points: will be increased to the nearest power of 2.
 template <typename FloatType>
-void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int in_outer_perimeter_points,
+void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int input_outer_perimeter_points,
                                                         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) {
 //  double Rho, Theta, Phi;
   if (from_Rho > to_Rho || from_Theta > to_Theta || from_Phi > to_Phi
-      || in_outer_perimeter_points < 1 || radius_points < 1)
+      || input_outer_perimeter_points < 1 || radius_points < 1)
     throw std::invalid_argument("Error! Invalid argument for RunFieldCalculationPolar() !");
-  int outer_perimeter_points = in_outer_perimeter_points;
-  if (outer_perimeter_points != 1) outer_perimeter_points = ceil_to_2_pow_n<FloatType>(in_outer_perimeter_points);
+  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 = LeRu_cutoff(std::complex<FloatType>(to_Rho, 0));
+  auto near_field_nmax = calcNmax(to_Rho);
+  SetMaxTerms(near_field_nmax);
 
   std::vector<std::vector<FloatType> >  Pi(outer_perimeter_points), Tau(outer_perimeter_points);
   for (auto &val:Pi) val.resize(near_field_nmax);
@@ -454,7 +454,6 @@ void MultiLayerMie<FloatType>::RunFieldCalculationPolar(const int in_outer_perim
   }
 
 
-
   // Calculate scattering coefficients an_ and bn_
 //  calcScattCoeffs();
   // Calculate expansion coefficients aln_,  bln_, cln_, and dln_

src/nmie.hpp

@@ -138,7 +138,7 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
     // Run calculation
     void RunMieCalculation();
     void RunFieldCalculation();
-    void RunFieldCalculationPolar(const int in_outer_perimeter_points = 1,
+    void RunFieldCalculationPolar(const int input_outer_perimeter_points = 1,
                                   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),
@@ -222,8 +222,8 @@ inline std::complex<T> my_exp(const std::complex<T>& x) {
     std::vector< std::vector<FloatType> > coords_;
 
   private:
-    void calcNstop();
-    void calcNmax(unsigned int first_layer);
+    unsigned int calcNstop(FloatType xL = -1);
+    unsigned int calcNmax(FloatType xL = -1);
 
     std::complex<FloatType> calc_an(int n, FloatType XL, std::complex<FloatType> Ha, std::complex<FloatType> mL,
                                  std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -44,7 +44,7 @@ std::vector< std::complex<double>>
 
 
 void parse_mpmath_data(const double min_abs_tol, const std::tuple< std::complex<double>, int, std::complex<double>, double, double > data,
-                       std::complex<double> &z, int &n, std::complex<double> &func_mp,
+                       std::complex<double> &z, unsigned int &n, std::complex<double> &func_mp,
                        double &re_abs_tol, double &im_abs_tol){
   z = std::get<0>(data);
   n = std::get<1>(data);
@@ -134,7 +134,7 @@ TEST(bn_test, mpmath_generated_input) {
 TEST(zeta_psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 2e-10;
   std::complex<double> z, zeta_mp;
-  int n;
+  unsigned int n;
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : zeta_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, zeta_mp, re_abs_tol, im_abs_tol);
@@ -196,7 +196,7 @@ TEST(zeta_psizeta_test, mpmath_generated_input) {
 TEST(psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 9e-11;
   std::complex<double> z, PsiZeta_mp;
-  int n;
+  unsigned int n;
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : psi_mul_zeta_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, PsiZeta_mp, re_abs_tol, im_abs_tol);
@@ -221,7 +221,7 @@ TEST(psizeta_test, mpmath_generated_input) {
 TEST(psi_test, mpmath_generated_input) {
   double min_abs_tol = 1e-12;
   std::complex<double> z, Psi_mp;
-  int n;
+  unsigned int n;
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : psi_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, Psi_mp, re_abs_tol, im_abs_tol);
@@ -243,7 +243,7 @@ TEST(psi_test, mpmath_generated_input) {
 TEST(D3test, mpmath_generated_input) {
   double min_abs_tol = 2e-11;
   std::complex<double> z, D3_mp;
-  int n;
+  unsigned int n;
   double re_abs_tol,  im_abs_tol;
   for (const auto &data : D3_test_16digits) {
     parse_mpmath_data(min_abs_tol, data, z, n, D3_mp, re_abs_tol, im_abs_tol);

tests/test_near_field.cc

@@ -15,7 +15,10 @@ TEST(RunFieldCalculationPolar, HandlesInput) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
   EXPECT_THROW(nmie.RunFieldCalculationPolar(0), std::invalid_argument);
   EXPECT_THROW(nmie.RunFieldCalculationPolar(1,1,10,5), std::invalid_argument);
-  nmie.RunFieldCalculationPolar();
+  nmie.SetLayersSize({0.099});
+  nmie.SetLayersIndex({ {0.75,0}});
+  nmie.RunMieCalculation();
+  nmie.RunFieldCalculationPolar(2);
 }
 //TEST(BulkSphere, HandlesInput) {
 //  nmie::MultiLayerMie<nmie::FloatType> nmie;