if d=0 then an matches for mesomie and nmie

src/mesomie.hpp

@@ -74,22 +74,37 @@ void MesoMie<FloatType>::calc_ab_classic(int nmax,
   evalPsiZetaD1D3(x * m,                       //
                   Psi_mx, Zeta_mx, D1_mx, D3_mx);
   for (int n = 0; n <= nmax; n++) {
+    // D1(x) = 0.573785843449477
+    // D1(mx) = (5.0024999985617e-7 - 0.99999999974975j)
+    // D3(x) = (-9.99900009998611e-7 + 0.999900009999j)
+    // psi_n = -0.867382528698782
+    // ksi_n = (-0.867382528698782 + 0.497742452386882j)
+    if (n < 3) {
+      std::cout << "n = " << n << std::endl;
+      std::cout << "D1(x) = " << D1_x[n] << std::endl;
+      std::cout << "D1(mx) = " << D1_mx[n] << std::endl;
+      std::cout << "D3(x) = " << D3_x[n] << std::endl;
+      std::cout << "psi(x) = " << Psi_x[n] << std::endl;
+      std::cout << "ksi(x) = " << Zeta_x[n] << std::endl;
+    }
     an_cl[n] = Psi_x[n] *
-               (                           //
-                   m * D1_x[n] - D1_mx[n]  //
-                   ) /                     //
-               Zeta_x[n] *
-               (                           //
-                   m * D3_x[n] - D1_mx[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]  //
-               );                          //
+               (                                           //
+                   D1_x[n] - m * D1_mx[n]                  //
+                   ) /                                     //
+               (                                           //
+                   Zeta_x[n] * (                           //
+                                   D3_x[n] - m * D1_mx[n]  //
+                                   )                       //
+               );                                          //
   }
 }
 
@@ -116,35 +131,38 @@ void MesoMie<FloatType>::calc_ab(int nmax,
 
   for (int n = 0; n <= nmax; n++) {
     an_[n] = Psi_xd[n] *
-             (                                                           //
-                 eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] +         //
-                 (eps_m - eps_d) *                                       //
-                     (                                                   //
-                         static_cast<FloatType>(n * (n + 1)) * d_perp +  //
-                         xd * D1_xd[n] * xm * D1_xm[n] * d_parallel      //
-                         ) /                                             //
-                     R                                                   //
-                 ) /                                                     //
-             Zeta_xd[n] *
-             (                                                           //
-                 eps_m * xd * D3_xd[n] - eps_d * xm * D1_xm[n] +         //
-                 (eps_m - eps_d) *                                       //
-                     (                                                   //
-                         static_cast<FloatType>(n * (n + 1)) * d_perp +  //
-                         xd * D3_xd[n] * xm * D1_xm[n] * d_parallel      //
-                         ) /                                             //
-                     R                                                   //
+             (                                                               //
+                 eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] +             //
+                 (eps_m - eps_d) *                                           //
+                     (                                                       //
+                         static_cast<FloatType>(n * (n + 1)) * d_perp +      //
+                         xd * D1_xd[n] * xm * D1_xm[n] * d_parallel          //
+                         ) /                                                 //
+                     R                                                       //
+                 ) /                                                         //
+             (                                                               //
+                 Zeta_xd[n] *                                                //
+                 (                                                           //
+                     eps_m * xd * D3_xd[n] - eps_d * xm * D1_xm[n] +         //
+                     (eps_m - eps_d) *                                       //
+                         (                                                   //
+                             static_cast<FloatType>(n * (n + 1)) * d_perp +  //
+                             xd * D3_xd[n] * xm * D1_xm[n] * d_parallel      //
+                             ) /                                             //
+                         R                                                   //
+                     )                                                       //
              );
     bn_[n] = Psi_xd[n] *
-             (                                         //
-                 xd * D1_xd[n] - xm * D1_xm[n] +       //
-                 (xm * xm - xd * xd) * d_parallel / R  //
-                 ) /                                   //
-             Zeta_xd[n] *
-             (                                         //
-                 xd * D3_xd[n] - xm * D1_xm[n] +       //
-                 (xm * xm - xd * xd) * d_parallel / R  //
-             );                                        //
+             (                                                          //
+                 xd * D1_xd[n] - xm * D1_xm[n] +                        //
+                 (xm * xm - xd * xd) * d_parallel / R                   //
+                 ) /                                                    //
+             (                                                          //
+                 Zeta_xd[n] * (                                         //
+                                  xd * D3_xd[n] - xm * D1_xm[n] +       //
+                                  (xm * xm - xd * xd) * d_parallel / R  //
+                                  )                                     //
+             );                                                         //
   }
 }
 

tests/mpmath_riccati_bessel.py

@@ -11,51 +11,88 @@ def LeRu_cutoff(z):
 # Wu, Wang, Radio Science, Volume 26, Number 6, Pages 1393-1401, November-December 1991,
 # after eq 13f.
 # Riccati-Bessel z*j_n(z)
-def psi(n,z):
-    return mp.sqrt( (mp.pi * z)/2 ) * mp.autoprec(mp.besselj)(n+1/2,z)
+def psi(n, z):
+    return mp.sqrt((mp.pi * z)/2) * mp.autoprec(mp.besselj)(n+1/2, z)
 # Riccati-Bessel -z*y_n(z)
-def xi(n,z):
-    return -mp.sqrt( (mp.pi * z)/2 ) * mp.autoprec(mp.bessely)(n+1/2,z)
+
+
+def xi(n, z):
+    return -mp.sqrt((mp.pi * z)/2) * mp.autoprec(mp.bessely)(n+1/2, z)
 # Riccati-Bessel psi - i* xi
-def ksi(n,z):
-    return psi(n,z) - 1.j * xi(n,z)
 
-def psi_div_ksi(n,z):
-    return psi(n,z)/ksi(n,z)
 
-def psi_mul_ksi(n,z):
-    return psi(n,z)*ksi(n,z)
+def ksi(n, z):
+    return psi(n, z) - 1.j * xi(n, z)
+
+
+def psi_div_ksi(n, z):
+    return psi(n, z)/ksi(n, z)
+
+
+def psi_mul_ksi(n, z):
+    return psi(n, z)*ksi(n, z)
+
 
-def psi_div_xi(n,z):
-    return psi(n,z)/xi(n,z)
+def psi_div_xi(n, z):
+    return psi(n, z)/xi(n, z)
 
 # to compare r(n,z) with Wolfram Alpha
 # n=49, z=1.3-2.1i,  SphericalBesselJ[n-1,z]/SphericalBesselJ[n,z]
-def r(n,z):
+
+
+def r(n, z):
     if n > 0:
-        return psi(n-1,z)/psi(n,z)
+        return psi(n-1, z)/psi(n, z)
     return mp.cos(z)/mp.sin(z)
 
 
 def D(n, z, f):
-    return f(n-1,z)/f(n,z) - n/z
+    return f(n-1, z)/f(n, z) - n/z
 
-def D1(n,z):
-    if n == 0: return mp.cos(z)/mp.sin(z)
+
+def D1(n, z):
+    if n == 0:
+        return mp.cos(z)/mp.sin(z)
     return D(n, z, psi)
 
 # Wolfram Alpha example D2(10, 10-10j): SphericalBesselY[9, 10-10i]/SphericalBesselY[10,10-10i]-10/(10-10i)
-def D2(n,z):
-    if n == 0: return -mp.sin(z)/mp.cos(z)
+
+
+def D2(n, z):
+    if n == 0:
+        return -mp.sin(z)/mp.cos(z)
     return D(n, z, xi)
 
-def D3(n,z):
-    if n == 0: return 1j
+
+def D3(n, z):
+    if n == 0:
+        return 1j
     return D(n, z, ksi)
 
 
 # bulk sphere
+
+# Le Ru
+def an_lr(n, x, m):
+    print(f'D1(x) = {D1(n, x)}\n\
+            D1(mx) = {D1(n, m*x)}\n\
+            D3(x) = {D3(n, x)}\n\
+            psi_n = {psi(n,x)}\n\
+            ksi_n = {ksi(n,x)}\n\
+    ')
+
+    return (
+        (
+            psi(n, x)*(m*D1(n, x)-D1(n, m*x))
+        ) /
+        (
+            ksi(n, x)*(m * D3(n, x)-D1(n, m*x))
+        )
+    )
+
 # Ovidio
+
+
 def an(n, x, m):
     # print(f'D1 = {D1(n, m*x)}\n\
     # psi_n = {psi(n,x)}\n\
@@ -64,9 +101,11 @@ def an(n, x, m):
     # 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) )
+        ((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))
     )
+
+
 def bn(n, x, m):
     # print(f'D1 = {D1(n, m*x)}\n\
     # psi_n = {psi(n,x)}\n\
@@ -75,8 +114,8 @@ def bn(n, x, m):
     # 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) )
+        ((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))
     )
 
 # # Du
@@ -93,34 +132,41 @@ def bn(n, x, m):
 #         (r(n,mx)*m*ksi(n,x)-ksi(n-1,x))
 #     )
 
-def Qext_diff(n,x,m):
-    return ((2*n +1)*2/x**2) * (an(n,x,m) + bn(n,x,m)).real
 
-def Qsca_diff(n,x,m):
-    return ((2*n +1)*2/x**2) * (abs(an(n,x,m))**2 + abs(bn(n,x,m))**2)
+def Qext_diff(n, x, m):
+    return ((2*n + 1)*2/x**2) * (an(n, x, m) + bn(n, x, m)).real
+
+
+def Qsca_diff(n, x, m):
+    return ((2*n + 1)*2/x**2) * (abs(an(n, x, m))**2 + abs(bn(n, x, m))**2)
+
 
 def Qany(x, m, nmax, output_dps, Qdiff):
     Qany = 0
     Qlist = []
-    Qprev=""
+    Qprev = ""
     convergence_max_length = 35
     i = 0
     for n in range(1, nmax+1):
-        diff = Qdiff(n,x,m)
+        diff = Qdiff(n, x, m)
         Qany += diff
         Qnext = mp.nstr(Qany, output_dps)
         Qlist.append(Qany)
         i += 1
-        if Qprev !=Qnext: i = 0
+        if Qprev != Qnext:
+            i = 0
         Qprev = Qnext
         print(n, ' ', end='', flush=True)
-        if i >= convergence_max_length: break
+        if i >= convergence_max_length:
+            break
     print('')
     Qlist.append(Qany)
     return Qlist
 
+
 def Qsca(x, m, nmax, output_dps):
     return Qany(x, m, nmax, output_dps, Qsca_diff)
 
+
 def Qext(x, m, nmax, output_dps):
     return Qany(x, m, nmax, output_dps, Qext_diff)

tests/test_bulk_sphere.cc

@@ -43,7 +43,7 @@ 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.099, {1.75, 0}, 7.417859e-06, 7.417859e-06, 'z'},
+        // {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'},
@@ -52,7 +52,7 @@ std::vector<std::tuple<double, std::complex<double>, double, double, char> >
         // {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'},
+        {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'},
@@ -107,7 +107,7 @@ TEST(BulkSphere, MesoMie) {
     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++) {
+    for (int n = 0; n < nmax && n < 200; n++) {
       std::cout << n << an_nmie[n] << ' ' << an_meso[n + 1] << ' '
                 << an_meso_cl[n + 1] << std::endl;
     }