scatt/tests/test_Riccati_Bessel_logarithmic_derivative.cc

#include "gtest/gtest.h"
#include "../src/nmie-impl.hpp"
#include "test_spec_functions_data.hpp"
// From W. Yang APPLIED OPTICS Vol. 42, No. 9,  20 March 2003
// Dtest refractive index is m={1.05,1}, the size parameter is x = 80
std::vector<int> Dtest_n({0,1,30,50,60,70,75,80,85,90,99,116,130});
std::vector< std::complex<double>>
    Dtest_D1({
      //Orig
//                 {0.11449e-15 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 },
//                 {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935},
//                 {0.13645,-0.10019e+01 },{0.18439,-0.10070e+01 },
//                 {0.21070,-0.10107e+01 },{0.23845,-0.10154e+01 },
//                 {0.26752,-0.10210e+01 },{0.29777,-0.10278e+01 },
//                 {0.35481,-0.10426e+01 },{0.46923,-0.10806e+01 },
//                 {0.17656,-0.13895e+01 }
      // mod (from Python mpmath)
                 {0.0,-1.0}, {7.464603828e-5,-0.9999958865},
                 {0.03476380918,-0.9986960672},{0.09529213152,-0.999347654},
                 {0.1364513887,-1.001895883},{0.184388335,-1.006979164},
                 {0.2107044267,-1.01072099},{0.2384524295,-1.015382914},
                 {0.2675164524,-1.021040337},{0.2977711192,-1.027753418},
                 {0.3548096904,-1.042622957},{0.4692294405,-1.080629479},
                 {0.5673827836,-1.121108944},
             });

std::vector< std::complex<double>>
    Dtest_D2({{0.64966e-69 ,-0.10000e+01 },{0.74646e-04 ,-0.10000e+01 },
              {0.34764e-01 ,-0.99870},{0.95292e-01 ,-0.99935},
              {0.13645,-0.10019e+01 },{0.17769,-0.10099e+01 },
              {0.41264e-01 ,-0.21076e+01 },{-0.20190,0.10435e+01 },
              {-0.26343,0.10223e+01 },{-0.29339,0.10291e+01 },
              {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 },
              {-0.56047,0.11206e+01 }});
std::vector< std::complex<double>>
    Dtest_D3({{0.00000,0.10000e+01 },{-0.73809e-04 ,0.10000e+01 },
              {-0.34344e-01 ,0.99912},{-0.94022e-01 ,0.10004e+01 },
              {-0.13455,0.10032e+01 },{-0.18172,0.10084e+01 },
              {-0.20762,0.10122e+01 },{-0.23494,0.10169e+01 },
              {-0.26357,0.10225e+01 },{-0.29339,0.10291e+01 },
              {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 },
              {-0.56047,0.11206e+01 }});


int LeRu_cutoff(std::complex<double> z) {
  auto x = std::abs(z);
  return std::round(x + 11 * std::pow(x, (1.0 / 3.0)) + 1);
}

void parse_mpmath_data(const double min_abs_tol, const std::tuple< std::complex<double>, int, std::complex<double>, double, double > data,
                       std::complex<double> &z, int &n, std::complex<double> &func_mp,
                       double &re_abs_tol, double &im_abs_tol){
  z = std::get<0>(data);
  n = std::get<1>(data);
  func_mp = std::get<2>(data);
  re_abs_tol = ( std::get<3>(data) > min_abs_tol && std::real(func_mp) < min_abs_tol)
                    ? std::get<3>(data) : min_abs_tol;
  im_abs_tol = ( std::get<4>(data) > min_abs_tol && std::imag(func_mp) < min_abs_tol)
                    ? std::get<4>(data) : min_abs_tol;
  // if re(func_mp) < 0.5 then round will give 0. To avoid zero tolerance add one.
  re_abs_tol *= std::abs(std::round(std::real(func_mp))) + 1;
  im_abs_tol *= std::abs(std::round(std::imag(func_mp))) + 1;
}
void parse2_mpmath_data(const nmie::FloatType min_abs_tol,
                        const std::tuple< nmie::FloatType, std::complex<nmie::FloatType>, int, std::complex<nmie::FloatType>, nmie::FloatType, nmie::FloatType > data,
                        nmie::FloatType &x, std::complex<nmie::FloatType> &m, int &n, std::complex<nmie::FloatType> &func_mp,
                        nmie::FloatType &re_abs_tol, nmie::FloatType &im_abs_tol){
  x = std::get<0>(data);
  m = std::get<1>(data);
  n = std::get<2>(data);
  func_mp = std::get<3>(data);
  re_abs_tol = ( std::get<4>(data) > min_abs_tol && std::real(func_mp) < min_abs_tol)
               ? std::get<4>(data) : min_abs_tol;
  im_abs_tol = ( std::get<5>(data) > min_abs_tol && std::imag(func_mp) < min_abs_tol)
               ? std::get<5>(data) : min_abs_tol;
  // if re(func_mp) < 0.5 then round will give 0. To avoid zero tolerance add one.
  re_abs_tol *= std::abs(std::round(std::real(func_mp))) + 1;
  im_abs_tol *= std::abs(std::round(std::imag(func_mp))) + 1;
}

template<class T> inline T pow2(const T value) {return value*value;}

//TEST(an_test, DISABLED_mpmath_generated_input) {
TEST(an_test, mpmath_generated_input) {
  double min_abs_tol = 3e-14, x;
  std::complex<double> m, an_mp;
  int n;
  double re_abs_tol,  im_abs_tol;
  for (const auto &data : an_test_30digits) {
    parse2_mpmath_data(min_abs_tol, data, x, m, n, an_mp, re_abs_tol, im_abs_tol);
    auto Nstop = LeRu_cutoff(m*x)+1;

    nmie::MultiLayerMie<nmie::FloatType> ml_mie;
    ml_mie.SetLayersSize({x});
    ml_mie.SetLayersIndex({m});
    ml_mie.SetMaxTerms(Nstop);
    ml_mie.calcScattCoeffs();
    auto an = ml_mie.GetAn();
//    auto bn = ml_mie.GetBn();

    if (n > an.size()) continue;
    if (n == 0) continue;
    EXPECT_NEAR(std::real(an[n-1]), std::real(an_mp), re_abs_tol)
              << "Db at n=" << n << " Nstop="<< Nstop<<" m="<<m<<" x="<<x;
    EXPECT_NEAR(std::imag(an[n-1]), std::imag(an_mp), im_abs_tol)
              << "Db at n=" << n << " Nstop="<< Nstop<<" m="<<m<<" x="<<x;
  }
}

TEST(zeta_psizeta_test, DISABLED_mpmath_generated_input) {
//TEST(zeta_psizeta_test, mpmath_generated_input) {
  double min_abs_tol = 2e-10;
  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+135), D3(Nstop+135),
        PsiZeta(Nstop+135), Psi(Nstop);
    nmie::evalDownwardD1(z, D1dr);
    nmie::evalUpwardD3(z, D1dr, D3, PsiZeta);
    nmie::evalUpwardPsi(z, D1dr, Psi);
    auto a = std::real(PsiZeta[n]);
    auto b = std::imag(PsiZeta[n]);
    auto c = std::real(Psi[n]);
    auto d = std::imag(Psi[n]);
    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::vector<std::complex<nmie::FloatType>> D1dr(Nstop+35), D3(Nstop+35), zeta(Nstop);
//    nmie::evalDownwardD1(z, D1dr);
//    nmie::evalUpwardD3(z, D1dr, D3);
//    nmie::evalUpwardZeta(z, D3, zeta);

    EXPECT_NEAR(std::real(zeta), std::real(zeta_mp), re_abs_tol)
              << "zeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
    EXPECT_NEAR(std::imag(zeta), std::imag(zeta_mp), im_abs_tol)
              << "zeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
  }
}

// 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, DISABLED_mpmath_generated_input) {
//TEST(psizeta_test, mpmath_generated_input) {
  double min_abs_tol = 9e-11;
  std::complex<double> z, PsiZeta_mp;
  int n;
  double re_abs_tol,  im_abs_tol;
  for (const auto &data : psi_mul_zeta_test_16digits) {
    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), D3(Nstop), PsiZeta(Nstop);
    nmie::evalDownwardD1(z, D1dr);
    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;
    EXPECT_NEAR(std::imag(PsiZeta[n]), std::imag(PsiZeta_mp), im_abs_tol)
              << "PsiZeta at n=" << n << " Nstop="<< Nstop<<" z="<<z;
//    std::vector<nmie::FloatType> PsiUp(Nstop);
//    nmie::evalPsi(std::real(z), PsiUp);
//    EXPECT_NEAR(((PsiUp[n])), std::real(PsiZeta_mp), re_abs_tol)
//              << "PsiZeta(up) at n=" << n << " z="<<z;
  }
}


TEST(psi_test, mpmath_generated_input) {
  double min_abs_tol = 1e-12;
  std::complex<double> z, Psi_mp;
  int n;
  double re_abs_tol,  im_abs_tol;
  for (const auto &data : psi_test_16digits) {
    parse_mpmath_data(min_abs_tol, data, z, n, Psi_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+35), Psi(Nstop);
    nmie::evalDownwardD1(z, D1dr);
    nmie::evalUpwardPsi(z, D1dr, Psi);

    EXPECT_NEAR(std::real(Psi[n]), std::real(Psi_mp), re_abs_tol)
              << "Psi at n=" << n << " Nstop="<< Nstop<<" z="<<z;
    EXPECT_NEAR(std::imag(Psi[n]), std::imag(Psi_mp), im_abs_tol)
              << "Psi at n=" << n << " Nstop="<< Nstop<<" z="<<z;
  }
}


TEST(D3test, DISABLED_mpmath_generated_input) {
//TEST(D3test, mpmath_generated_input) {
  double min_abs_tol = 2e-11;
  std::complex<double> z, D3_mp;
  int n;
  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)+35;
    std::vector<std::complex<nmie::FloatType>> D1dr(Nstop), D3(Nstop), PsiZeta(Nstop);
    nmie::evalDownwardD1(z, D1dr);
    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;
    EXPECT_NEAR(std::imag(D3[n]), std::imag(D3_mp), im_abs_tol)
              << "D3 at n=" << n << " Nstop="<< Nstop<<" z="<<z;

  }
}


//TEST(D1test, DISABLED_mpmath_generated_input) {
  TEST(D1test, mpmath_generated_input) {
  double min_abs_tol = 2e-11, x;
  std::complex<double> m, z, D1_mp;
  int n;
  double re_abs_tol,  im_abs_tol;
  for (const auto &data : D1_test_30digits) {
    parse2_mpmath_data(min_abs_tol, data, x, m, n, D1_mp, re_abs_tol, im_abs_tol);
    z = m*x;
    auto Nstop = LeRu_cutoff(z)+1;
    std::vector<std::complex<nmie::FloatType>> Db(Nstop),Dold(Nstop+135), r;
    int valid_digits = 14;
    int nstar = nmie::getNStar(Nstop, z, valid_digits);
    r.resize(nstar);
    nmie::evalBackwardR(z,r);
    nmie::convertRtoD1(z, r, Db);
    if (n > Db.size()) continue;
    EXPECT_NEAR(std::real(Db[n]), std::real(D1_mp), re_abs_tol)
              << "Db at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
    EXPECT_NEAR(std::imag(Db[n]), std::imag(D1_mp), im_abs_tol)
              << "Db at n=" << n << " Nstop="<< Nstop<<" nstar="<<nstar<< " z="<<z;
    nmie::evalDownwardD1(z, Dold);
    if (n > Dold.size()) continue;
    EXPECT_NEAR(std::real(Dold[n]), std::real(D1_mp), re_abs_tol)
              << "Dold at n=" << n << " Nstop="<< Nstop<< " z="<<z;
    EXPECT_NEAR(std::imag(Dold[n]), std::imag(D1_mp), im_abs_tol)
              << "Dold at n=" << n << " Nstop="<< Nstop<< " z="<<z;

  }
}


//TEST(D1test, DISABLED_WYang_data){
TEST(D1test, WYang_data){
  double abs_tol = 1e-9;
  int test_loss_digits = std::round(15 - std::log10(1/abs_tol));
  int Nstop = 131;
  std::vector<std::complex<nmie::FloatType>> Df(Nstop), Db(Nstop),Dold(Nstop), r;
  std::complex<nmie::FloatType> z(1.05,1);
  z = z*80.0;
// eval D1 directly from backward recurrence
  nmie::evalDownwardD1(z, Dold);

//  eval forward recurrence
  r.resize(Nstop+1);
  nmie::evalForwardR(z, r);
  nmie::convertRtoD1(z, r, Df);
// eval backward recurrence
  int valid_digits = 6;
  int nstar = nmie::getNStar(Nstop, z, valid_digits);
  r.resize(nstar);
  nmie::evalBackwardR(z,r);
  nmie::convertRtoD1(z, r, Db);

for (int i = 0; i < Dtest_n.size(); i++) {
  int n = Dtest_n[i];
  int forward_loss_digits = nmie::evalKapteynNumberOfLostSignificantDigits(n, z);
  forward_loss_digits += 3; // Kapteyn is too optimistic
  if (test_loss_digits > forward_loss_digits ) {
    EXPECT_NEAR(std::real(Df[n]), std::real(Dtest_D1[i]),
                abs_tol) << "f at n=" << n << " lost digits = " << forward_loss_digits;
    EXPECT_NEAR(std::imag(Df[n]), std::imag(Dtest_D1[i]),
                abs_tol) << "f at n=" << n << " lost digits = " << forward_loss_digits;
  }
  EXPECT_NEAR(std::real(Db[n]), std::real(Dtest_D1[i]),
              abs_tol) << "b at n=" << n;
  EXPECT_NEAR(std::imag(Db[n]), std::imag(Dtest_D1[i]),
              abs_tol) << "b at n=" << n;
  if (n < Dold.size()-15) {
    EXPECT_NEAR(std::real(Dold[n]), std::real(Dtest_D1[i]),
                abs_tol) << "old at n=" << n;
    EXPECT_NEAR(std::imag(Dold[n]), std::imag(Dtest_D1[i]),
                abs_tol) << "old at n=" << n;
  }

}
}



TEST(KaptyenTest, HandlesInput) {
  // H.Du APPLIED OPTICS, Vol. 43, No. 9, 20 March 2004
  double l = nmie::evalKapteynNumberOfLostSignificantDigits(80, std::complex<double>(100,100));
  EXPECT_EQ(l, 7)<<"Should be equal";
  std::complex<double> z(10000,0);
  l = nmie::evalKapteynNumberOfLostSignificantDigits(5070, z);
  EXPECT_EQ(l, 0)<<"Should be equal";
  // find NStar such that l_nstar(z) - l_nmax(z) >= valid_digits
  int NStar = nmie::getNStar(5070, z,6);
  EXPECT_GE(NStar, 10130);
//  const double pi=3.14159265358979323846;
//  z = std::complex<double>(100,100);
//  l = nmie::evalKapteynNumberOfLostSignificantDigits(1, z);
//  EXPECT_EQ(l, 0)<<"Should be equal";

}


int main(int argc, char **argv) {
  testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}