Du test passes for mesomie

src/mesomie.hpp

@@ -59,65 +59,37 @@
 namespace nmie {
 //******************************************************************************
 template <typename FloatType>
-void MesoMie<FloatType>::calc_ab_classic(int nmax,
-                                         FloatType x,
-                                         std::complex<FloatType> m) {
-  an_cl.resize(nmax + 1, static_cast<FloatType>(0.0));
-  bn_cl.resize(nmax + 1, static_cast<FloatType>(0.0));
-  std::vector<std::complex<FloatType>>         //
-      Psi_mx(nmax + 1), Zeta_mx(nmax + 1),     //
-      Psi_x(nmax + 1), Zeta_x(nmax + 1),       //
-      D1_mx(nmax + 1), D3_mx(nmax + 1),        //
-      D1_x(nmax + 1), D3_x(nmax + 1);          //
-  evalPsiZetaD1D3(std::complex<FloatType>(x),  //
-                  Psi_x, Zeta_x, D1_x, D3_x);  //
-  evalPsiZetaD1D3(x * m,                       //
-                  Psi_mx, Zeta_mx, D1_mx, D3_mx);
-  for (int n = 0; n <= nmax; n++) {
-    // D1(x) = 0.573785843449477
-    // D1(mx) = (5.0024999985617e-7 - 0.99999999974975j)
-    // D3(x) = (-9.99900009998611e-7 + 0.999900009999j)
-    // psi_n = -0.867382528698782
-    // ksi_n = (-0.867382528698782 + 0.497742452386882j)
-    if (n < 3) {
-      std::cout << "n = " << n << std::endl;
-      std::cout << "D1(x) = " << D1_x[n] << std::endl;
-      std::cout << "D1(mx) = " << D1_mx[n] << std::endl;
-      std::cout << "D3(x) = " << D3_x[n] << std::endl;
-      std::cout << "psi(x) = " << Psi_x[n] << std::endl;
-      std::cout << "ksi(x) = " << Zeta_x[n] << std::endl;
-    }
-    an_cl[n] = Psi_x[n] *
-               (                                           //
-                   m * D1_x[n] - D1_mx[n]                  //
-                   ) /                                     //
-               (                                           //
-                   Zeta_x[n] * (                           //
-                                   m * D3_x[n] - D1_mx[n]  //
-                                   )                       //
-               );
-    bn_cl[n] = Psi_x[n] *
-               (                                           //
-                   D1_x[n] - m * D1_mx[n]                  //
-                   ) /                                     //
-               (                                           //
-                   Zeta_x[n] * (                           //
-                                   D3_x[n] - m * D1_mx[n]  //
-                                   )                       //
-               );                                          //
+void MesoMie<FloatType>::calc_Q() {
+  auto nmax = an_.size();
+  Qext_ = 0.0;
+  Qsca_ = 0.0;
+
+  for (int n = nmax - 2; n >= 1; n--) {
+    //      for (int n = 0; n < nmax_; n++) {
+    const int n1 = n;
+    // Equation (27)
+    Qext_ += (n1 + n1 + 1.0) * (an_[n].real() + bn_[n].real());
+    // std::cout << n1 << ": " << Qext_ << "  ";
+    // Equation (28)
+    Qsca_ += (n1 + n1 + 1.0) *
+             (an_[n].real() * an_[n].real() + an_[n].imag() * an_[n].imag() +
+              bn_[n].real() * bn_[n].real() + bn_[n].imag() * bn_[n].imag());
   }
+  Qext_ *= 2 / pow2(x_);
+  Qsca_ *= 2 / pow2(x_);
 }
 
 //******************************************************************************
 template <typename FloatType>
-void MesoMie<FloatType>::calc_ab(int nmax,
-                                 FloatType R,
+void MesoMie<FloatType>::calc_ab(FloatType R,
                                  std::complex<FloatType> xd,
                                  std::complex<FloatType> xm,
                                  std::complex<FloatType> eps_d,
                                  std::complex<FloatType> eps_m,
                                  std::complex<FloatType> d_parallel,
                                  std::complex<FloatType> d_perp) {
+  x_ = R;
+  int nmax = std::round(x_ + 11 * std::pow(x_, (1.0 / 3.0)) + 1);
   an_.resize(nmax + 1, static_cast<FloatType>(0.0));
   bn_.resize(nmax + 1, static_cast<FloatType>(0.0));
   std::vector<std::complex<FloatType>>      //

src/nmie.hpp

@@ -534,15 +534,15 @@ template <typename FloatType = double>
 class MesoMie {
  public:
   std::vector<std::complex<FloatType>> an_, bn_;
+  FloatType x_;
+  std::complex<FloatType> m_;
   std::vector<std::complex<FloatType>> GetAn() { return an_; };
   std::vector<std::complex<FloatType>> GetBn() { return bn_; };
 
-  std::vector<std::complex<FloatType>> an_cl, bn_cl;
-  std::vector<std::complex<FloatType>> GetAnClassic() { return an_cl; };
-  std::vector<std::complex<FloatType>> GetBnClassic() { return bn_cl; };
-
-  void calc_ab(int nmax,
-               FloatType R,
+  FloatType Qsca_ = 0.0, Qext_ = 0.0;
+  FloatType GetQsca() { return Qsca_; };
+  FloatType GetQext() { return Qext_; };
+  void calc_ab(FloatType R,
                std::complex<FloatType> xd,
                std::complex<FloatType> xm,
                std::complex<FloatType> eps_d,
@@ -550,7 +550,12 @@ class MesoMie {
                std::complex<FloatType> d_parallel,
                std::complex<FloatType> d_perp);
 
-  void calc_ab_classic(int nmax, FloatType x, std::complex<FloatType> m);
+  void calc_Q();
+  // template <typename outputType = FloatType>
+  // outputType GetQext();
+
+  // template <typename outputType = FloatType>
+  // outputType GetQsca();
 };  // end of class MesoMie
 
 }  // end of namespace nmie

tests/test_bulk_sphere.cc

@@ -42,20 +42,19 @@ TEST(BulkSphere, ArgPi) {
 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.099, {1.75, 0}, 7.417859e-06, 7.417859e-06, 'z'},
-        // {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'},
+        {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'},
+        {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'},
     };
 //******************************************************************************
 TEST(BulkSphere, DISABLED_HandlesInput) {
@@ -88,39 +87,25 @@ TEST(BulkSphere, MesoMie) {
   for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
     auto m = std::get<1>(data);
-    //    auto Nstop = nmie::LeRu_near_field_cutoff(m*x)+1;
-
-    nmie.SetLayersSize({x});
-    nmie.SetLayersIndex({m});
-    //    nmie.SetMaxTerms(Nstop);
-    nmie.calcScattCoeffs();
-    auto an_nmie = nmie.GetAn();
-    int nmax = an_nmie.size();
-    mesomie.calc_ab(nmax + 2,       // nmax
-                    x,              // R
-                    {x, 0},         // xd
-                    x * m,          // xm
-                    1.0 / (x * m),  // eps_d
-                    m / x,          // eps_m
-                    {0, 0},         // d_parallel
-                    {0, 0});        // d_perp
-    auto an_meso = mesomie.GetAn();
-    mesomie.calc_ab_classic(nmax + 2, x, m);
-    auto an_meso_cl = mesomie.GetAnClassic();
-    for (int n = 0; n < nmax && n < 200; n++) {
-      std::cout << n << an_nmie[n] << ' ' << an_meso[n + 1] << ' '
-                << an_meso_cl[n + 1] << std::endl;
-    }
+    mesomie.calc_ab(x,        // R
+                    {x, 0},   // xd
+                    x * m,    // xm
+                    {1, 0},   // eps_d
+                    m * m,    // eps_m
+                    {0, 0},   // d_parallel
+                    {0, 0});  // d_perp
+                              // eps_m * xd / (eps_d * xm)
+    mesomie.calc_Q();
 
-    // 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);
+    double Qext = static_cast<double>(mesomie.GetQext());
+    double Qsca = static_cast<double>(mesomie.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) {