Zeta is now evaluated from eq 18

src/nmie-impl.hpp

@@ -392,8 +392,9 @@ namespace nmie {
   void MultiLayerMie<FloatType>::calcD1D3(const std::complex<FloatType> z,
                                std::vector<std::complex<FloatType> >& D1,
                                std::vector<std::complex<FloatType> >& D3) {
+    std::vector<std::complex<FloatType> > PsiZeta(nmax_+1);
     evalDownwardD1(z, D1);
-    evalUpwardD3 (z, D1, D3);
+    evalUpwardD3 (z, D1, D3, PsiZeta);
   }
 
 
@@ -413,12 +414,18 @@ namespace nmie {
   void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
                                   std::vector<std::complex<FloatType> >& Psi,
                                   std::vector<std::complex<FloatType> >& Zeta) {
-    std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1);
+    std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1),
+        PsiZeta(nmax_+1);
     // First, calculate the logarithmic derivatives
-    calcD1D3(z, D1, D3);
-    // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
+    evalDownwardD1(z, D1);
+    // Now, use the upward recurrence to calculate Psi equations (20ab)
     evalUpwardPsi(z,  D1, Psi);
-    evalUpwardZeta(z, D3, Zeta);
+    // Now, use the upward recurrence to calculate Psi*Zeta equations (18ad)
+    evalUpwardD3 (z, D1, D3, PsiZeta);
+    for (int i = 0; i < Zeta.size(); i++) {
+      Zeta[i] = PsiZeta[i]/Psi[i];
+    }
+//    evalUpwardZeta(z, D3, Zeta);
   }
 
 

src/special-functions-impl.hpp

@@ -337,28 +337,9 @@ void evalDownwardD1 (const std::complex<FloatType> z,
 
 void evalUpwardD3 (const std::complex<FloatType> z,
                    const std::vector<std::complex<FloatType> >& D1,
-                   std::vector<std::complex<FloatType> >& D3) {
-  int nmax = D1.size()-1;
-  std::vector<std::complex<FloatType> > PsiZeta(nmax+1);
-
-  // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
-  PsiZeta[0] = static_cast<FloatType>(0.5)*(static_cast<FloatType>(1.0) - std::complex<FloatType>(nmm::cos(2.0*z.real()), nmm::sin(2.0*z.real()))
-      *static_cast<FloatType>(nmm::exp(-2.0*z.imag())));
-  D3[0] = std::complex<FloatType>(0.0, 1.0);
-  const std::complex<FloatType> zinv = std::complex<FloatType>(1.0, 0.0)/z;
-  for (int n = 1; n <= nmax; n++) {
-    PsiZeta[n] = PsiZeta[n - 1]*(static_cast<FloatType>(n)*zinv - D1[n - 1])
-        *(static_cast<FloatType>(n)*zinv - D3[n - 1]);
-    D3[n] = D1[n] + std::complex<FloatType>(0.0, 1.0)/PsiZeta[n];
-  }
-}
-
-void evalUpwardPsiZeta (const std::complex<FloatType> z,
-                   const std::vector<std::complex<FloatType> >& D1,
+                   std::vector<std::complex<FloatType> >& D3,
                    std::vector<std::complex<FloatType> >& PsiZeta) {
   int nmax = D1.size()-1;
-  std::vector<std::complex<FloatType> > D3(nmax+1);
-
   // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
   PsiZeta[0] = static_cast<FloatType>(0.5)*(static_cast<FloatType>(1.0) - std::complex<FloatType>(nmm::cos(2.0*z.real()), nmm::sin(2.0*z.real()))
       *static_cast<FloatType>(nmm::exp(-2.0*z.imag())));
@@ -371,6 +352,7 @@ void evalUpwardPsiZeta (const std::complex<FloatType> z,
   }
 }
 
+
 void evalUpwardPsi (const std::complex<FloatType> z,
                     const std::vector<std::complex<FloatType> > D1,
                    std::vector<std::complex<FloatType> >& Psi) {
@@ -396,36 +378,36 @@ void evalUpwardZeta (const std::complex<FloatType> z,
   }
 }
 
-void evalForwardRiccatiBessel(const FloatType x, const FloatType first, const FloatType second,
-                              std::vector<FloatType> &values) {
-  values[0] = first;
-  values[1] = second;
-  int nmax = values.size();
-  for (int i = 1; i < nmax-1; i++) {
-    values[i+1] = (1 + 2*i) * values[i]/x - values[i-1];
-  }
-}
-
-void evalChi(const FloatType x, std::vector<FloatType> &Chi) {
-  auto first = nmm::cos(x);
-  auto second = first/x + nmm::sin(x);
-  evalForwardRiccatiBessel(x, first, second, Chi);
-}
-
-void evalPsi(const FloatType x, std::vector<FloatType> &Psi) {
-  auto first = nmm::sin(x);
-  auto second = first/x - nmm::cos(x);
-  evalForwardRiccatiBessel(x, first, second, Psi);
-}
-
-void composeZeta(const std::vector<FloatType> &Psi,
-                 const std::vector<FloatType> &Chi,
-                 std::vector< std::complex<FloatType>> &Zeta) {
-  int nmax = Zeta.size();
-  for (int i = 0; i < nmax; i++) {
-    Zeta[i] = std::complex<FloatType > (Psi[i], Chi[i]);
-  }
-}
+//void evalForwardRiccatiBessel(const FloatType x, const FloatType first, const FloatType second,
+//                              std::vector<FloatType> &values) {
+//  values[0] = first;
+//  values[1] = second;
+//  int nmax = values.size();
+//  for (int i = 1; i < nmax-1; i++) {
+//    values[i+1] = (1 + 2*i) * values[i]/x - values[i-1];
+//  }
+//}
+//
+//void evalChi(const FloatType x, std::vector<FloatType> &Chi) {
+//  auto first = nmm::cos(x);
+//  auto second = first/x + nmm::sin(x);
+//  evalForwardRiccatiBessel(x, first, second, Chi);
+//}
+//
+//void evalPsi(const FloatType x, std::vector<FloatType> &Psi) {
+//  auto first = nmm::sin(x);
+//  auto second = first/x - nmm::cos(x);
+//  evalForwardRiccatiBessel(x, first, second, Psi);
+//}
+//
+//void composeZeta(const std::vector<FloatType> &Psi,
+//                 const std::vector<FloatType> &Chi,
+//                 std::vector< std::complex<FloatType>> &Zeta) {
+//  int nmax = Zeta.size();
+//  for (int i = 0; i < nmax; i++) {
+//    Zeta[i] = std::complex<FloatType > (Psi[i], Chi[i]);
+//  }
+//}
 
 }  // end of namespace nmie
 #endif  // SRC_SPECIAL_FUNCTIONS_IMPL_HPP_

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -64,8 +64,8 @@ void parse_mpmath_data(const double min_abs_tol, const std::tuple< std::complex<
 
 template<class T> inline T pow2(const T value) {return value*value;}
 
-//TEST(zeta_test, DISABLED_mpmath_generated_input) {
-TEST(zeta_test, mpmath_generated_input) {
+//TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
+TEST(zeta_psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 2e-5;
   std::complex<double> z, zeta_mp;
   int n;
@@ -74,10 +74,10 @@ TEST(zeta_test, mpmath_generated_input) {
     parse_mpmath_data(min_abs_tol, data, z, n, zeta_mp, re_abs_tol, im_abs_tol);
     auto Nstop = LeRu_cutoff(z)+10000;
     if (n > Nstop) continue;
-    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+135),
+    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+135), D3(Nstop+135),
         PsiZeta(Nstop+135), Psi(Nstop);
     nmie::evalDownwardD1(z, D1dr);
-    nmie::evalUpwardPsiZeta(z, D1dr, PsiZeta);
+    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
     nmie::evalUpwardPsi(z, D1dr, Psi);
     auto a = std::real(PsiZeta[n]);
     auto b = std::imag(PsiZeta[n]);
@@ -86,8 +86,8 @@ TEST(zeta_test, mpmath_generated_input) {
     auto c_one = std::complex<nmie::FloatType>(0, 1);
     auto zeta = (a*c + b*d)/(pow2(c) + pow2(d)) +
         c_one * ((b*c - a*d)/(pow2(c) + pow2(d)));
+//    zeta = PsiZeta[n]/Psi[n];
     if (std::isnan(std::real(zeta)) || std::isnan(std::imag(zeta))) continue;
-//    std::complex<nmie::FloatType > zeta = PsiZeta[n]/Psi[n];
 //    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+35), D3(Nstop+35), zeta(Nstop);
 //    nmie::evalDownwardD1(z, D1dr);
 //    nmie::evalUpwardD3(z, D1dr, D3);
@@ -100,6 +100,31 @@ TEST(zeta_test, mpmath_generated_input) {
   }
 }
 
+// Old way to evaluate Zeta
+TEST(zeta_test, DISABLED_mpmath_generated_input) {
+//TEST(zeta_test, mpmath_generated_input) {
+  double min_abs_tol = 2e-5;
+  std::complex<double> z, zeta_mp;
+  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);
+    auto Nstop = LeRu_cutoff(z)+10000;
+    if (n > Nstop) continue;
+    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop),
+        PsiZeta(Nstop), zeta(Nstop);
+    nmie::evalDownwardD1(z, D1dr);
+    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
+    nmie::evalUpwardZeta(z, D3, zeta);
+    if (std::isnan(std::real(zeta[n])) || std::isnan(std::imag(zeta[n]))) continue;
+
+    EXPECT_NEAR(std::real(zeta[n]), std::real(zeta_mp), re_abs_tol)
+              << "zeta[n] at n=" << n << " Nstop="<< Nstop<<" z="<<z;
+    EXPECT_NEAR(std::imag(zeta[n]), std::imag(zeta_mp), im_abs_tol)
+              << "zeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
+  }
+}
+
 
 TEST(psizeta_test, mpmath_generated_input) {
   double min_abs_tol = 9e-11;
@@ -110,9 +135,9 @@ TEST(psizeta_test, mpmath_generated_input) {
     parse_mpmath_data(min_abs_tol, data, z, n, PsiZeta_mp, re_abs_tol, im_abs_tol);
     auto Nstop = LeRu_cutoff(z)+10000;
     if (n > Nstop) continue;
-    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+135), PsiZeta(Nstop+135);
+    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
     nmie::evalDownwardD1(z, D1dr);
-    nmie::evalUpwardPsiZeta(z, D1dr, PsiZeta);
+    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
 
     EXPECT_NEAR(std::real(PsiZeta[n]), std::real(PsiZeta_mp), re_abs_tol)
               << "PsiZeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
@@ -155,10 +180,10 @@ TEST(D3test, mpmath_generated_input) {
   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);
-    auto Nstop = LeRu_cutoff(z)+1;
-    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop+35), D3(Nstop+35);
+    auto Nstop = LeRu_cutoff(z)+35;
+    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
     nmie::evalDownwardD1(z, D1dr);
-    nmie::evalUpwardD3(z, D1dr, D3);
+    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
 
     EXPECT_NEAR(std::real(D3[n]), std::real(D3_mp), re_abs_tol)
               << "D3 at n=" << n << " Nstop="<< Nstop<<" z="<<z;

tests/test_bulk_sphere.cc

@@ -57,10 +57,12 @@ TEST(BulkSphere, HandlesInput) {
     nmie.RunMieCalculation();
     double Qext = static_cast<double>(nmie.GetQext());
     double Qsca = static_cast<double>(nmie.GetQsca());
-    EXPECT_FLOAT_EQ(std::get<2>(data), Qext)
+    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(std::get<3>(data), Qsca)
+    EXPECT_FLOAT_EQ(Qsca_Du, Qsca)
               << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
   }