an test passes

src/nmie.cc

@@ -74,7 +74,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m,
+  int ScattCoeffs(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m,
                   const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
 
     if (x.size() != L || m.size() != L)
@@ -118,7 +118,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
+  int ExpanCoeffs(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax,
                   std::vector<std::vector<std::complex<double> > >& aln, std::vector<std::vector<std::complex<double> > >& bln,
                   std::vector<std::vector<std::complex<double> > >& cln, std::vector<std::vector<std::complex<double> > >& dln) {
 

src/nmie.hpp

@@ -42,11 +42,19 @@
 #include <boost/math/constants/constants.hpp>
 #endif
 namespace nmie {
-  int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn);
+  int ScattCoeffs(const unsigned int L, const int pl,
+                  const std::vector<double>& x, const std::vector<std::complex<double> >& m,
+                  const int nmax,
+                  std::vector<std::complex<double> >& an,
+                  std::vector<std::complex<double> >& bn);
 
-  int ExpanCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax,
-                std::vector<std::vector<std::complex<double> > >& an, std::vector<std::vector<std::complex<double> > >& bn,
-                std::vector<std::vector<std::complex<double> > >& cn, std::vector<std::vector<std::complex<double> > >& dn);
+  int ExpanCoeffs(const unsigned int L, const int pl,
+                  const std::vector<double>& x, const std::vector<std::complex<double> >& m,
+                  const int nmax,
+                  std::vector<std::vector<std::complex<double> > >& an,
+                  std::vector<std::vector<std::complex<double> > >& bn,
+                  std::vector<std::vector<std::complex<double> > >& cn,
+                  std::vector<std::vector<std::complex<double> > >& dn);
 
 //helper functions
 

tests/mpmath_riccati_bessel.py

@@ -57,6 +57,12 @@ def D3(n,z):
 # bulk sphere
 # Ovidio
 def an(n, x, m):
+    print(f'D1 = {D1(n, m*x)}\n\
+    psi_n = {psi(n,x)}\n\
+    psi_nm1 = {psi(n-1,x)}\n\
+    ksi_n = {ksi(n,x)}\n\
+    ksi_nm1 = {ksi(n-1,x)}\n\
+    ')
     return (
             ( ( D1(n, m*x)/m + n/x )*psi(n,x) - psi(n-1,x) ) /
             ( ( D1(n, m*x)/m + n/x )*ksi(n,x) - ksi(n-1,x) )

tests/mpmath_special_functions_test_generator.py

@@ -111,7 +111,10 @@ class UpdateSpecialFunctionsEvaluations:
             mpf_m = mp.mpc(mr, mi)
             z = mpf_x*mpf_m
             if self.is_only_x: z = mp.mpf(x)
-            mpf_value = func(n, z)
+            if self.is_xm:
+                mpf_value = func(n, mpf_x, mpf_m)
+            else:
+                mpf_value = func(n, z)
             z_str = self.compose_result_string(mpf_x, mpf_m, n, mpf_value, output_dps)
             if mp.nstr(mpf_value.real, output_dps) == '0.0' \
                     or mp.nstr(mpf_value.imag, output_dps) == '0.0':
@@ -159,8 +162,9 @@ class UpdateSpecialFunctionsEvaluations:
         output_list.append('};')
         return output_list
 
-    def run_test(self, func, funcname, is_only_x=False):
+    def run_test(self, func, funcname, is_only_x=False, is_xm=False):
         self.is_only_x = is_only_x
+        self.is_xm = is_xm
         self.remove_argument_duplicates()
         out_list_result = self.get_test_data(self.complex_arguments, self.output_dps,
                                              self.max_num_elements_of_nlist,
@@ -181,7 +185,7 @@ class UpdateSpecialFunctionsEvaluations:
                     filtered_list[-1][1] == self.complex_arguments[i][1]):
                 continue
             # argument list is sorted, so when only x is needed
-            # the record with the largest m
+            # keep the record with the largest m
             if (self.is_only_x
                     and filtered_list[-1][0] == self.complex_arguments[i][0]):
                 # continue
@@ -197,15 +201,18 @@ def main():
     sf_evals = UpdateSpecialFunctionsEvaluations(filename='test_spec_functions_data.hpp',
                                                  complex_arguments=mia.complex_arguments,
                                                  output_dps=30, max_num_elements_of_nlist=51)
+                                                 # output_dps=7, max_num_elements_of_nlist=51)
                                                  # output_dps=5, max_num_elements_of_nlist=3)
-    sf_evals.run_test(mrb.D1, 'D1')
+    # sf_evals.run_test(mrb.D1, 'D1')
     # sf_evals.run_test(mrb.D2, 'D2')
     # sf_evals.run_test(mrb.D3, 'D3')
     # sf_evals.run_test(mrb.psi, 'psi', is_only_x=True)
-    # # In literature Zeta or Ksi denote the Riccati-Bessel function of third kind.
     # sf_evals.run_test(mrb.xi, 'xi', is_only_x=True)
+    # # In literature Zeta or Ksi denote the Riccati-Bessel function of third kind.
     # sf_evals.run_test(mrb.ksi, 'zeta', is_only_x=True)
 
+    # sf_evals.run_test(mrb.an, 'an', is_xm=True)
+
     # sf_evals.run_test(mrb.psi, 'psi')
     # sf_evals.run_test(mrb.psi_div_ksi, 'psi_div_ksi')
     # sf_evals.run_test(mrb.psi_mul_ksi, 'psi_mul_zeta', is_only_x=True)
@@ -213,23 +220,25 @@ def main():
     with open(sf_evals.filename, 'w') as out_file:
         out_file.write(sf_evals.get_file_content())
 
-    # for record in mia.complex_arguments:
-    #     mp.mp.dps = 20
-    #     output_dps = 7
-    #     x = mp.mpf(str(record[0]))
-    #     mr = str(record[1][0])
-    #     mi = str(record[1][1])
-    #     m = mp.mpc(mr, mi)
-    #     Qext_ref = record[2]
-    #     Qsca_ref = record[3]
-    #     test_case = record[4]
-    #     nmax = int(x + 4.05*x**(1./3.) + 2)+2+28
-    #     print(f"\n ===== test case: {test_case} =====", flush=True)
-    #     print(f"x={x}, m={m}, N={nmax} \nQsca_ref = {Qsca_ref}    \tQext_ref = {Qext_ref}", flush=True)
-    #     Qext_mp = mrb.Qext(x,m,nmax, output_dps)
-    #     Qsca_mp = mrb.Qsca(x,m,nmax, output_dps)
-    #     print(f"Qsca_mp  = {mp.nstr(Qsca_mp[-1],output_dps)}    \tQext_mp  = {mp.nstr(Qext_mp[-1],output_dps)}", flush=True)
-    #     print(mp.nstr(Qsca_mp,output_dps))
-    #     print(mp.nstr(Qext_mp,output_dps))
+    for record in mia.complex_arguments:
+        mp.mp.dps = 20
+        output_dps = 16
+        x = mp.mpf(str(record[0]))
+        mr = str(record[1][0])
+        mi = str(record[1][1])
+        m = mp.mpc(mr, mi)
+        Qext_ref = record[2]
+        Qsca_ref = record[3]
+        test_case = record[4]
+        nmax = int(x + 4.05*x**(1./3.) + 2)+2+28
+        print(f"\n ===== test case: {test_case} =====", flush=True)
+        print(f"x={x}, m={m}, N={nmax} \nQsca_ref = {Qsca_ref}    \tQext_ref = {Qext_ref}", flush=True)
+        # Qext_mp = mrb.Qext(x,m,nmax, output_dps)
+        # Qsca_mp = mrb.Qsca(x,m,nmax, output_dps)
+        # print(f"Qsca_mp  = {mp.nstr(Qsca_mp[-1],output_dps)}    \tQext_mp  = {mp.nstr(Qext_mp[-1],output_dps)}", flush=True)
+        # print(mp.nstr(Qsca_mp,output_dps))
+        # print(mp.nstr(Qext_mp,output_dps))
+        # n=1
+        # print(f'n={n}, x={x}, m={m}\nan[{n}]={mp.nstr(mrb.an(n,x,m), output_dps)}')
 
 main()

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -80,6 +80,33 @@ void parse2_mpmath_data(const nmie::FloatType min_abs_tol,
 
 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;

tests/test_spec_functions_data.hpp

@@ -1652,6 +1652,50 @@ D3_test_16digits
 {{3.141592653589793,0.0},20,{-6.285084985730819,5.832688628116996e-28},6.3e-16,5.8e-44},
 };
 
+// x, complex(m), n, complex(f(n,z)), abs_err_real, abs_err_imag
+std::vector< std::tuple< nmie::FloatType, std::complex<nmie::FloatType>, int, std::complex<nmie::FloatType>, nmie::FloatType, nmie::FloatType > >
+an_test_30digits
+= {
+{100.0,{1.33,0.00001},0,{0.983292861475793547449358879534,0.124840935725346993952308423791},9.8e-31,1.2e-31},
+{100.0,{1.33,0.00001},1,{0.984860446612096957241775578166,-0.117441925244364000775108414221},9.8e-31,1.2e-31},
+{100.0,{1.33,0.00001},2,{0.985965322504954123969621952786,0.114030053006644283488612943716},9.9e-31,1.1e-31},
+{100.0,{1.33,0.00001},3,{0.98309422036283488576271509011,-0.124415977608877776986379398593},9.8e-31,1.2e-31},
+{100.0,{1.33,0.00001},4,{0.991371042866121608625494977157,0.08795829271464134280082195607},9.9e-31,8.8e-32},
+{100.0,{1.33,0.00001},5,{0.980248634667341669556147870197,-0.134834855608139093882796207474},9.8e-31,1.3e-31},
+{100.0,{1.33,0.00001},6,{0.99718055122375359392774331865,0.0449047618462165643766318875834},1.0e-30,4.5e-32},
+{100.0,{1.33,0.00001},7,{0.977109178727649399333768793089,-0.145371593321445950228768159659},9.8e-31,1.5e-31},
+{100.0,{1.33,0.00001},8,{0.998944269454026267720243481287,-0.0168214926226047272238686782359},1.0e-30,1.7e-32},
+{100.0,{1.33,0.00001},9,{0.97483843031762829521227929905,-0.15244923102990775776867544793},9.7e-31,1.5e-31},
+{100.0,{1.33,0.00001},10,{0.989682373029210168246874877077,-0.0972464738249357923933102153541},9.9e-31,9.7e-32},
+{100.0,{1.33,0.00001},11,{0.974058276029222943006235123659,-0.154763341782793180850316831588},9.7e-31,1.5e-31},
+{100.0,{1.33,0.00001},13,{0.97321520791555822224253070484,-0.157411119794584445279087903581},9.7e-31,1.6e-31},
+{100.0,{1.33,0.00001},14,{0.904663248735925755303113778303,-0.292361654042369775061262128286},9.0e-31,2.9e-31},
+{100.0,{1.33,0.00001},16,{0.821789077184583765787653734649,-0.381603517691708934777072956747},8.2e-31,3.8e-31},
+{100.0,{1.33,0.00001},18,{0.72613260329953407155848417394,-0.444886204316627242871792444943},7.3e-31,4.4e-31},
+{100.0,{1.33,0.00001},20,{0.646163397269584115804313346767,-0.47700354344730323479394821332},6.5e-31,4.8e-31},
+{100.0,{1.33,0.00001},22,{0.608259288093308282827444683986,-0.486845852056768417880388842449},6.1e-31,4.9e-31},
+{100.0,{1.33,0.00001},24,{0.597202161138536973717709609856,-0.489173914624679657860238877562},6.0e-31,4.9e-31},
+{100.0,{1.33,0.00001},27,{0.230834302937511370779284381301,-0.420262916810272826582474221898},2.3e-31,4.2e-31},
+{100.0,{1.33,0.00001},30,{0.183770761359558719881501512678,-0.386322169474363897617474153971},1.8e-31,3.9e-31},
+{100.0,{1.33,0.00001},34,{0.00407268980589458002534477977901,0.0556762984125619636521173223988},4.1e-33,5.6e-32},
+{100.0,{1.33,0.00001},38,{0.0472729617907451310636812562269,0.209550902880260799159827077153},4.7e-32,2.1e-31},
+{100.0,{1.33,0.00001},42,{0.466131617991045553592382970539,0.49810420892156931844106101529},4.7e-31,5.0e-31},
+{100.0,{1.33,0.00001},47,{0.966899506343116582867783282697,0.176777123382178820318858700892},9.7e-31,1.8e-31},
+{100.0,{1.33,0.00001},52,{0.718207294717742029675948752342,-0.448958051273150730138669586336},7.2e-31,4.5e-31},
+{100.0,{1.33,0.00001},58,{0.0236414565408137004489218151802,-0.149452411551900295640136802398},2.4e-32,1.5e-31},
+{100.0,{1.33,0.00001},65,{0.790181580208819746691144800806,0.40597144972645670117974018671},7.9e-31,4.1e-31},
+{100.0,{1.33,0.00001},72,{0.189913560500779200402138612099,-0.391136586515346529172705738098},1.9e-31,3.9e-31},
+{100.0,{1.33,0.00001},80,{0.999005778595894389161376387959,-0.0141888952885494395739093887755},1.0e-30,1.4e-32},
+{100.0,{1.33,0.00001},89,{0.906221854741476619292691065599,0.290379310943678733900129380023},9.1e-31,2.9e-31},
+{100.0,{1.33,0.00001},99,{0.00660940947386932922358492734156,0.0775093664967703347987311255936},6.6e-33,7.8e-32},
+{100.0,{1.33,0.00001},111,{0.00000321366061193370200254328420947,0.00103983157486592840043188559885},3.2e-36,1.0e-33},
+{100.0,{1.33,0.00001},123,{5.17583806523184440915572503511e-13,2.67274951899386243525029046053e-10},5.2e-43,2.7e-40},
+{100.0,{1.33,0.00001},137,{4.71748803281812297397012823952e-24,-1.04259994939657750111183308868e-19},4.7e-54,1.0e-49},
+{100.0,{1.33,0.00001},153,{4.04619812001120256031833116666e-37,-1.26394242328878445327736126568e-32},4.0e-67,1.3e-62},
+{100.0,{1.33,0.00001},170,{6.54854383559658349678357032904e-53,-2.27878735404540029090834523813e-48},6.5e-83,2.3e-78},
+{100.0,{1.33,0.00001},190,{9.68143517313324082875587812679e-74,-3.55321407196406464604668533115e-69},9.7e-104,3.6e-99},
+};
+
 // complex(z), n, complex(f(n,z)), abs_err_real, abs_err_imag
 std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >
 psi_mul_zeta_test_16digits