Start of using bessel::

bessel.cc

@@ -105,7 +105,7 @@ namespace nmie {
       csy.resize(n+1);
       cdy.resize(n+1);
       std::complex<double> cf, cf0, cf1, cs, csa, csb;
-      int k, m;
+      int m;
       a0 = std::abs(z);
       nm = n;      
       if (a0 < 1.0e-60) {
@@ -209,7 +209,7 @@ namespace nmie {
       // !
     int msta2 ( double x, int n, int mp ) {
       double a0, ejn, f, f0, f1, hmp;
-      int it, n0, n1, nn;
+      int  n0, n1, nn;
       double obj;
       a0 = std::abs ( x );
       hmp = 0.5 * mp;
@@ -281,7 +281,7 @@ namespace nmie {
 // !
     int  msta1 ( double x, int mp ) {
       double a0, f, f0, f1;
-      int it, n0, n1, nn;
+      int n0, n1, nn;
       a0 = std::abs ( x );
       n0 = static_cast<int> (1.1 * a0 ) + 1;
       f0 = envj ( n0, a0 ) - mp;

compare.cc

@@ -250,18 +250,18 @@ int main(int argc, char *argv[]) {
     int samples = range.size();
     std::vector<double> zero(samples, 0.0);
     std::vector<double> Xp, Yp, Zp;
-    // // X line
-    // Xp.insert(Xp.end(), range.begin(), range.end());
-    // Yp.insert(Yp.end(), zero.begin(), zero.end());
-    // Zp.insert(Zp.end(), zero.begin(), zero.end());
+    // X line
+    Xp.insert(Xp.end(), range.begin(), range.end());
+    Yp.insert(Yp.end(), zero.begin(), zero.end());
+    Zp.insert(Zp.end(), zero.begin(), zero.end());
     // Y line
     Xp.insert(Xp.end(), zero.begin(), zero.end());
     Yp.insert(Yp.end(), range.begin(), range.end());
     Zp.insert(Zp.end(), zero.begin(), zero.end());
-    // // // Z line
-    // Xp.insert(Xp.end(), zero.begin(), zero.end());
-    // Yp.insert(Yp.end(), zero.begin(), zero.end());
-    // Zp.insert(Zp.end(), range.begin(), range.end());
+    // Z line
+    Xp.insert(Xp.end(), zero.begin(), zero.end());
+    Yp.insert(Yp.end(), zero.begin(), zero.end());
+    Zp.insert(Zp.end(), range.begin(), range.end());
     int ncoord = Xp.size();
     // Test solid sphere
     x = {size};

go.sh

@@ -15,7 +15,7 @@ echo Uncomment next line to compile compare.cc
 #  g++ -Ofast -std=c++11 compare.cc nmie.cc bessel.cc nmie-wrapper.cc -lm -lrt -o scattnlay-g.bin -ltcmalloc -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -g
 
 #DEBUG!
-clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 compare.cc nmie.cc bessel.cc nmie-old.cc -lm -lrt -o scattnlay.bin
+#clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 compare.cc nmie.cc bessel.cc nmie-old.cc -lm -lrt -o scattnlay.bin
 
 cp scattnlay.bin ../scattnlay
 # cp scattnlay-g.bin ../scattnlay-g
@@ -28,7 +28,7 @@ cp scattnlay.bin ../scattnlay
 # done
 
 PROGRAM='./scattnlay.bin'
-time  ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4  $PROGRAM -l 5 0.4642 1.8000 1.7000 0.7114 0.8000 0.7000 0.7393 1.2000 0.0900 0.9168 2.8000 0.2000 1.0000 1.5000 0.4000  -t 0.0 90.0 5 -c test01
+#time  ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4  $PROGRAM -l 5 0.4642 1.8000 1.7000 0.7114 0.8000 0.7000 0.7393 1.2000 0.0900 0.9168 2.8000 0.2000 1.0000 1.5000 0.4000  -t 0.0 90.0 5 -c test01
 # #result
 # # test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01
 
@@ -94,11 +94,11 @@ time  ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4  $PROGRAM -l 5 0.4642 1.
 # ./$file.bin
 
 # ### Cython
-# rm scattnlay.so
-# export CFLAGS='-std=c++11'
-# python setup.py build_ext --inplace
-# cp scattnlay.so tests/python/
-# cd tests/python/
-# ./field.py
+rm scattnlay.so
+export CFLAGS='-std=c++11'
+python setup.py build_ext --inplace
+cp scattnlay.so tests/python/
+cd tests/python/
+./field.py
 # ./test01.py
 # ./test01-wrapper.py

nmie.cc

@@ -649,6 +649,12 @@ namespace nmie {
                    *(static_cast<double>(n)*zinv- D3[n - 1]);
       D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
     }
+    std::vector<std::complex<double> >  ZetaZ(nmax_+1), dZetaZ(nmax_ + 1);
+    bessel::calcZeta(ZetaZ, dZetaZ, nmax_, z);
+    for (int n = 0; n < nmax_+1; ++n) {
+      D3[n]=dZetaZ[n]/ZetaZ[n];
+    }
+
   }
 
 
@@ -681,6 +687,10 @@ namespace nmie {
       Psi[n] = Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
       Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
     }
+    
+    std::vector<std::complex<double> >  dZetaZ(nmax_ + 1);
+    bessel::calcZeta(Zeta, dZetaZ, nmax_, z);
+
 
   }
 
@@ -1188,10 +1198,12 @@ namespace nmie {
 
     // Check the result and change  aln_[0][n] and aln_[0][n] for exact zero
     for (int n = 0; n < nmax_; ++n) {
-      if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
-      else throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
-      if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
+      printf("n=%d, aln_=%g,%g,   bln_=%g,%g \n", n, real(aln_[0][n]), imag(aln_[0][n]),
+	     real(bln_[0][n]), imag(bln_[0][n]));
+      if (std::abs(aln_[0][n]) < 1e-1) aln_[0][n] = 0.0;
       else throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
+      if (std::abs(bln_[0][n]) < 1e-1) bln_[0][n] = 0.0;
+      else throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
     }
 
     isIntCoeffsCalc_ = true;

setup.py

@@ -53,7 +53,7 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
       license = 'GPL',
       platforms = 'any',
       ext_modules = [Extension("scattnlay",
-                               ["nmie.cc", "py_nmie.cc", "scattnlay.cpp"],
+                               ["nmie.cc", "bessel.cc","py_nmie.cc", "scattnlay.cpp"],
                                language = "c++",
                                include_dirs = [np.get_include()])], 
       extra_compile_args=['-std=c++11']

setup_cython.py

@@ -54,7 +54,7 @@ O. Pena, U. Pal, Comput. Phys. Commun. 180 (2009) 2348-2354.""",
       license = 'GPL',
       platforms = 'any',
       ext_modules = cythonize("scattnlay.pyx",                                        # our Cython source
-                              sources = ["nmie.cc", "py_nmie.cc", "scattnlay.cpp"],   # additional source file(s)
+                              sources = ["nmie.cc", "bessel.cc", "py_nmie.cc", "scattnlay.cpp"],   # additional source file(s)
                               language = "c++",                                       # generate C++ code
                               extra_compile_args = ['-std=c++11'],
       )

tests/python/field-dielectric-sphere.py

@@ -66,8 +66,8 @@ terms, E, H = fieldnlay(x, m, coord)
 Er = np.absolute(E)
 
 # |E|/|Eo|
-#Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
-Eh = Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2
+Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+#Eh = Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2
 
 result = np.vstack((coordX, coordY, coordZ, Eh)).transpose()
 
@@ -91,13 +91,13 @@ try:
     min_tick = max(0.1, min(min_tick, np.amin(edata)))
     max_tick = max(max_tick, np.amax(edata))
     #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
-    scale_ticks = np.linspace(min_tick,min_tick*1.2, 11)
+    scale_ticks = np.linspace(min_tick,max_tick, 11)
     #scale_ticks = np.linspace(0, 2, 11)
 
     # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
     cax = ax.imshow(edata, interpolation = 'nearest', cmap = cm.afmhot,
-                    #origin = 'lower', vmin = min_tick, vmax = max_tick,
-                    origin = 'lower', vmin = 0.25, vmax = 1,
+                    origin = 'lower', vmin = min_tick, vmax = max_tick,
+                    #origin = 'lower', vmin = 0.25, vmax = 1,
                     extent = (min(scale_x), max(scale_x), min(scale_y), max(scale_y))
                     #,norm = LogNorm()
                     )