added equation for bulk sphere from Le Ru papaer

src/mesomie.hpp

@@ -57,12 +57,48 @@
 #include "special-functions-impl.hpp"
 
 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++) {
+    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]  //
+               );                          //
+  }
+}
+
 //******************************************************************************
+template <typename FloatType>
 void MesoMie<FloatType>::calc_ab(int nmax,
                                  FloatType R,
-                                 FloatType xd,
-                                 FloatType xm,
+                                 std::complex<FloatType> xd,
+                                 std::complex<FloatType> xm,
                                  std::complex<FloatType> eps_d,
                                  std::complex<FloatType> eps_m,
                                  std::complex<FloatType> d_parallel,
@@ -78,7 +114,7 @@ void MesoMie<FloatType>::calc_ab(int nmax,
   evalPsiZetaD1D3(std::complex<FloatType>(xd), Psi_xd, Zeta_xd, D1_xd, D3_xd);
   evalPsiZetaD1D3(std::complex<FloatType>(xm), Psi_xm, Zeta_xm, D1_xm, D3_xm);
 
-  for (int n = 1; n <= nmax; n++) {
+  for (int n = 0; n <= nmax; n++) {
     an_[n] = Psi_xd[n] *
              (                                                           //
                  eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] +         //

src/nmie.hpp

@@ -537,14 +537,20 @@ class MesoMie {
   std::vector<std::complex<FloatType>> GetAn() { return an_; };
   std::vector<std::complex<FloatType>> GetBn() { return bn_; };
 
-  void calc_ab(int n,
+  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 xd,
-               FloatType xm,
+               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);
+
+  void calc_ab_classic(int nmax, FloatType x, std::complex<FloatType> m);
 };  // end of class MesoMie
 
 }  // end of namespace nmie

tests/test_bulk_sphere.cc

@@ -42,23 +42,24 @@ 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.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'},
+        // {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'},
+        // {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'},
     };
 //******************************************************************************
-// TEST(BulkSphere, DISABLED_HandlesInput) {
-TEST(BulkSphere, HandlesInput) {
+TEST(BulkSphere, DISABLED_HandlesInput) {
+  // TEST(BulkSphere, HandlesInput) {
   nmie::MultiLayerMie<nmie::FloatType> nmie;
   for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
@@ -95,11 +96,20 @@ TEST(BulkSphere, MesoMie) {
     nmie.calcScattCoeffs();
     auto an_nmie = nmie.GetAn();
     int nmax = an_nmie.size();
-    mesomie.calc_ab(nmax + 2, 1, 2, x, {1, 0}, std::sqrt(m), {0, 0}, {0, 0});
+    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++) {
-      if (an_meso[n + 1].real() > 0.1)
-        std::cout << an_nmie[n] << ' ' << an_meso[n + 1] << std::endl;
+      std::cout << n << an_nmie[n] << ' ' << an_meso[n + 1] << ' '
+                << an_meso_cl[n + 1] << std::endl;
     }
 
     // double Qext = static_cast<double>(nmie.GetQext());