initial test (with runtime error)

src/nmie.hpp

@@ -534,6 +534,8 @@ template <typename FloatType = double>
 class MesoMie {
  public:
   std::vector<std::complex<FloatType>> an_, bn_;
+  std::vector<std::complex<FloatType>> GetAn() { return an_; };
+  std::vector<std::complex<FloatType>> GetBn() { return bn_; };
 
   void calc_ab(int n,
                FloatType R,

src/special-functions-impl.hpp

@@ -426,28 +426,6 @@ void evalUpwardD3(const std::complex<FloatType> z,
 
 //******************************************************************************
 template <typename FloatType>
-void evalPsiZetaD1D3(const std::complex<FloatType> cxd,
-                     std::vector<std::complex<FloatType>>& Psi,
-                     std::vector<std::complex<FloatType>>& Zeta,
-                     std::vector<std::complex<FloatType>>& D1,
-                     std::vector<std::complex<FloatType>>& D3) {
-  int nmax = Psi.size() - 1;
-  std::vector<std::complex<FloatType>> PsiZeta(nmax + 1);
-
-  for (int n = 0; n <= nmax; n++) {
-    D1[n] = std::complex<FloatType>(0.0, -1.0);
-    D3[n] = std::complex<FloatType>(0.0, 1.0);
-  }
-
-  evalDownwardD1(cxd, D1);
-  evalUpwardPsi(cxd, D1, Psi);
-  evalUpwardD3(cxd, D1, D3, PsiZeta);
-  for (unsigned int i = 0; i < Zeta.size(); i++) {
-    Zeta[i] = PsiZeta[i] / Psi[i];
-  }
-}
-//******************************************************************************
-template <typename FloatType>
 std::complex<FloatType> complex_sin(const std::complex<FloatType> z) {
   auto a = z.real();
   auto b = z.imag();
@@ -526,6 +504,28 @@ void evalUpwardZeta(const std::complex<FloatType> z,
 //     Zeta[i] = std::complex<FloatType > (Psi[i], Chi[i]);
 //   }
 // }
+//******************************************************************************
+template <typename FloatType>
+void evalPsiZetaD1D3(const std::complex<FloatType> cxd,
+                     std::vector<std::complex<FloatType>>& Psi,
+                     std::vector<std::complex<FloatType>>& Zeta,
+                     std::vector<std::complex<FloatType>>& D1,
+                     std::vector<std::complex<FloatType>>& D3) {
+  int nmax = Psi.size() - 1;
+  std::vector<std::complex<FloatType>> PsiZeta(nmax + 1);
+
+  for (int n = 0; n <= nmax; n++) {
+    D1[n] = std::complex<FloatType>(0.0, -1.0);
+    D3[n] = std::complex<FloatType>(0.0, 1.0);
+  }
+
+  evalDownwardD1(cxd, D1);
+  evalUpwardPsi(cxd, D1, Psi);
+  evalUpwardD3(cxd, D1, D3, PsiZeta);
+  for (unsigned int i = 0; i < Zeta.size(); i++) {
+    Zeta[i] = PsiZeta[i] / Psi[i];
+  }
+}
 
 }  // end of namespace nmie
 #endif  // SRC_SPECIAL_FUNCTIONS_IMPL_HPP_

tests/test_bulk_sphere.cc

@@ -1,5 +1,6 @@
-#include "../src/nmie-basic.hpp"
+#include <complex>
 #include "../src/mesomie.hpp"
+#include "../src/nmie-basic.hpp"
 #include "gtest/gtest.h"
 
 // TODO fails for MP with 100 digits. And 16 digits, which should be equal to
@@ -34,42 +35,35 @@ TEST(BulkSphere, ArgPi) {
 #endif
 
 //******************************************************************************
+// 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'},
+    };
+//******************************************************************************
 // TEST(BulkSphere, DISABLED_HandlesInput) {
 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_near_field_cutoff(m*x)+1;
-
     nmie.SetLayersSize({x});
     nmie.SetLayersIndex({m});
     //    nmie.SetMaxTerms(Nstop);
@@ -89,34 +83,7 @@ TEST(BulkSphere, HandlesInput) {
 //******************************************************************************
 TEST(BulkSphere, MesoMie) {
   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'},
-      };
+  nmie::MesoMie<nmie::FloatType> mesomie;
   for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
     auto m = std::get<1>(data);
@@ -125,16 +92,25 @@ TEST(BulkSphere, MesoMie) {
     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);
+    nmie.calcScattCoeffs();
+    auto an_nmie = nmie.GetAn();
+    int nmax = an_nmie.size();
+    mesomie.calc_ab(nmax + 2, 1, 1, x, 1, std::sqrt(m).real(), 0, 0);
+    // auto an_meso = mesomie.GetAn();
+    for (int n = 0; n < nmax; n++) {
+      if (an_nmie[n].real() > 0.1)
+        std::cout << an_nmie[n] << std::endl;
+    }
+
+    // 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) {