turned all special funciton definitions into templates + misc

src/nmie-basic.hpp

@@ -274,7 +274,9 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
 
-  int LeRu_cutoff(std::complex<double> z) {
+  template <typename FloatType>
+  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);
 //    return 10000;
@@ -295,7 +297,7 @@ namespace nmie {
       nmax_ = newround(xL + 4.0*pow(xL, 1.0/3.0) + 2);
     }
     //Le Ru
-    auto Nstop = nmie::LeRu_cutoff(static_cast<double>(xL))+1;
+    auto Nstop = nmie::LeRu_cutoff(std::complex<FloatType>(xL,0))+1;
     if (Nstop > nmax_) nmax_ = Nstop;
   }
 

src/nmie-nearfield.hpp

@@ -405,5 +405,60 @@ namespace nmie {
       }
     }  // end of for all field coordinates
   }  //  end of MultiLayerMie::RunFieldCalculation()
+
+template <typename FloatType>
+int ceil_to_2_pow_n(const int input_n) {
+  int n = 2;
+  while (input_n > n) n *= 2;
+  return 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));
+  // 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!");
+  return delta;
+}
+
+
+// input parameters:
+//         in_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,
+                                                        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)
+    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);
+//  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));
+
+  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]);
+  }
+
+
+
+  // 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.cc

@@ -97,7 +97,7 @@ namespace nmie {
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
     }
-  }
+  }  //end of int ScattCoeffs(...)
 
   //**********************************************************************************//
   // This function emulates a C call to calculate the expansion coefficients          //
@@ -147,7 +147,7 @@ namespace nmie {
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
     }
-  }
+  }  // end of int ExpanCoeffs(...)
 
 
   //**********************************************************************************//
@@ -222,8 +222,7 @@ namespace nmie {
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
     }
-  }
-
+  }  // end of int nMie(...)
 
   //**********************************************************************************//
   // This function is just a wrapper to call the full 'nMie' function with fewer      //
@@ -256,9 +255,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,
+           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) {
     return nmie::nMie(L, pl, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2, -1, -1);
-
   }
   //**********************************************************************************//
   // This function is just a wrapper to call the full 'nMie' function with fewer      //
@@ -426,5 +428,5 @@ namespace nmie {
       std::cerr << "Invalid argument: " << ia.what() << std::endl;
       throw std::invalid_argument(ia);
     }
-  }
+  }  // end of int nField(...)
 }  // end of namespace nmie

src/nmie.hpp

@@ -138,6 +138,12 @@ 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,
+                                  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));
+
     void calcScattCoeffs();
     void calcExpanCoeffs();
 

src/special-functions-impl.hpp

@@ -56,11 +56,9 @@
 #include "nmie-precision.hpp"
 
 namespace nmie {
-////helper functions
-//template<class T> inline T pow2(const T value) {return value*value;}
-
 // Note, that Kapteyn seems to be too optimistic (at least by 3 digits
 // in some cases) for forward recurrence, see D1test with WYang_data
+template <typename FloatType>
 int evalKapteynNumberOfLostSignificantDigits(const int ni,
                                              const std::complex<FloatType> zz) {
   using std::abs;  using std::imag; using std::real; using std::log; using std::sqrt; using std::round;
@@ -74,6 +72,7 @@ int evalKapteynNumberOfLostSignificantDigits(const int ni,
       )))/ log(10.));
 }
 
+template <typename FloatType>
 int getNStar(int nmax, std::complex<FloatType> z, const int valid_digits) {
   if (nmax == 0) nmax = 1;
   int nstar = nmax;
@@ -92,8 +91,8 @@ int getNStar(int nmax, std::complex<FloatType> z, const int valid_digits) {
 
 // Custom implementation of complex cot function to avoid overflow
 // if Im(z) < 0, then it evaluates cot(z) as conj(cot(conj(z)))
-const std::complex<FloatType>
-complex_cot(const std::complex<FloatType> z) {
+template <typename FloatType>
+std::complex<FloatType> complex_cot(const std::complex<FloatType> z) {
   auto Remx = z.real();
   auto Immx = z.imag();
   int sign = (Immx>0) ? 1: -1; // use complex conj if needed for exp and return
@@ -109,6 +108,7 @@ 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) {
   if (r.size() < 1) throw std::invalid_argument(
@@ -123,6 +123,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) {
   if (r.size() < 1) throw std::invalid_argument(
@@ -137,6 +138,7 @@ void evalBackwardR (const std::complex<FloatType> z,
 //  nmm::cout << "R0 = " << r[0] <<" at arg = "<<z<<'\n';
 }
 
+template <typename FloatType>
 void convertRtoD1(const std::complex<FloatType> z,
                   std::vector<std::complex<FloatType> >& r,
                   std::vector<std::complex<FloatType> >& D1) {
@@ -151,6 +153,7 @@ void convertRtoD1(const std::complex<FloatType> z,
 }
 
 // ********************************************************************** //
+template <typename FloatType>
 void evalForwardD (const std::complex<FloatType> z,
                      std::vector<std::complex<FloatType> >& D) {
   int nmax = D.size();
@@ -162,6 +165,7 @@ void evalForwardD (const std::complex<FloatType> z,
 }
 
 // ********************************************************************** //
+template <typename FloatType>
 void evalForwardD1 (const std::complex<FloatType> z,
                    std::vector<std::complex<FloatType> >& D) {
   if (D.size()<1) throw std::invalid_argument("Should have a leas one element!\n");
@@ -319,6 +323,7 @@ void evalForwardD1 (const std::complex<FloatType> z,
 // Output parameters:                                                               //
 //   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions                //
 //**********************************************************************************//
+template <typename FloatType>
 void evalDownwardD1 (const std::complex<FloatType> z,
                      std::vector<std::complex<FloatType> >& D1) {
   int nmax = D1.size() - 1;
@@ -344,6 +349,7 @@ void evalDownwardD1 (const std::complex<FloatType> z,
 }
 
 
+template <typename FloatType>
 void evalUpwardD3 (const std::complex<FloatType> z,
                    const std::vector<std::complex<FloatType> >& D1,
                    std::vector<std::complex<FloatType> >& D3,
@@ -362,6 +368,7 @@ void evalUpwardD3 (const std::complex<FloatType> z,
 }
 
 
+template <typename FloatType>
 void evalUpwardPsi (const std::complex<FloatType> z,
                     const std::vector<std::complex<FloatType> > D1,
                    std::vector<std::complex<FloatType> >& Psi) {
@@ -374,6 +381,7 @@ void evalUpwardPsi (const std::complex<FloatType> z,
 }
 
 // Sometimes in literature Zeta is also denoted as Ksi, it is a Riccati-Bessel function of third kind.
+template <typename FloatType>
 void evalUpwardZeta (const std::complex<FloatType> z,
                     const std::vector<std::complex<FloatType> > D3,
                     std::vector<std::complex<FloatType> >& Zeta) {

tests/CMakeLists.txt

@@ -8,6 +8,13 @@ set(LIBS ${LIBS} pthread)
 
 # -- Output tests in directory
 
+add_executable("test_near_field"
+        test_near_field.cc)
+target_link_libraries("test_near_field" ${LIBS})
+add_test(NAME "test_near_field"
+        COMMAND "test_near_field")
+
+
 # In included file test_spec_functions_data.hpp there are results of multiple
 # precision computation that may overflow double precision at compile time.
 set_source_files_properties(test_Riccati_Bessel_logarithmic_derivative.cc

tests/test_near_field.cc

@@ -0,0 +1,67 @@
+#include "gtest/gtest.h"
+#include "../src/nmie-basic.hpp"
+#include "../src/nmie-nearfield.hpp"
+
+TEST(ceil_to_2_pow_n, HandlesInput) {
+  EXPECT_EQ(16, nmie::ceil_to_2_pow_n<nmie::FloatType>(10));
+  EXPECT_EQ(32, nmie::ceil_to_2_pow_n<nmie::FloatType>(20));
+  EXPECT_EQ(32, nmie::ceil_to_2_pow_n<nmie::FloatType>(30));
+  EXPECT_EQ(128, nmie::ceil_to_2_pow_n<nmie::FloatType>(100));
+  EXPECT_EQ(256, nmie::ceil_to_2_pow_n<nmie::FloatType>(200));
+  EXPECT_EQ(128, nmie::ceil_to_2_pow_n<nmie::FloatType>(128));
+}
+
+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();
+}
+//TEST(BulkSphere, HandlesInput) {
+//  nmie::MultiLayerMie<nmie::FloatType> nmie;
+//  // A list of tests for a bulk sphere from
+//  // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
+//  // table 1: sphere size and refractive index
+//  // followed by resulting extinction and scattering efficiencies
+//  std::vector< std::tuple< double, std::complex<double>, double, double, char > >
+//      parameters_and_results
+//      {
+//          // x, {Re(m), Im(m)}, Qext, Qsca, test_name
+//          {0.099, {0.75,0}, 7.417859e-06, 7.417859e-06, 'a'},
+//          {0.101, {0.75,0}, 8.033538e-06, 8.033538e-06, 'b'},
+//          {10,    {0.75,0},     2.232265, 2.232265, 'c'},
+//          {100,   {1.33,1e-5}, 2.101321, 2.096594, 'e'},
+//          {0.055, {1.5, 1},    0.10149104, 1.131687e-05, 'g'},
+//          {0.056, {1.5, 1},   0.1033467, 1.216311e-05, 'h'},
+//          {100,   {1.5, 1},    2.097502, 1.283697, 'i'},
+//          {1,     {10,  10},   2.532993, 2.049405, 'k'},
+//          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
+//          {100,   {10,  10,},  2.071124, 1.836785, 'l'},
+//          {10000, {1.33,1e-5}, 2.004089, 1.723857, 'f'},
+//          {10000, {1.5, 1},    2.004368, 1.236574, 'j'},
+//          {10000, {10,  10},   2.005914, 1.795393, 'm'},
+//      };
+//  for (const auto &data : parameters_and_results) {
+//    auto x = std::get<0>(data);
+//    auto m = std::get<1>(data);
+////    auto Nstop = nmie::LeRu_cutoff(m*x)+1;
+//
+//    nmie.SetLayersSize({x});
+//    nmie.SetLayersIndex({m});
+////    nmie.SetMaxTerms(Nstop);
+//    nmie.RunMieCalculation();
+//    double Qext = static_cast<double>(nmie.GetQext());
+//    double Qsca = static_cast<double>(nmie.GetQsca());
+//    double Qext_Du =  std::get<2>(data);
+//    double Qsca_Du =  std::get<3>(data);
+//    EXPECT_FLOAT_EQ(Qext_Du, Qext)
+//              << "Extinction of the bulk sphere, test case:" << std::get<4>(data)
+//              << "\nnmax_ = " << nmie.GetMaxTerms();
+//    EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
+//              << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
+//  }
+//}
+int main(int argc, char **argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}