correct calculations of zeta and diff(zeta)

Makefile

@@ -22,10 +22,10 @@ source:
 cython: scattnlay.pyx
 	cython --cplus scattnlay.pyx
 
-python_ext: nmie.cc py_nmie.cc scattnlay.cpp
+python_ext: nmie.cc bessel.cc py_nmie.cc scattnlay.cpp
 	export CFLAGS='-std=c++11' && python setup.py build_ext --inplace
 
-cython_ext: nmie.cc py_nmie.cc scattnlay.pyx
+cython_ext: nmie.cc bessel.cc py_nmie.cc scattnlay.pyx
 	export CFLAGS='-std=c++11' && python setup_cython.py build_ext --inplace
 
 install:
@@ -43,8 +43,8 @@ builddeb:
 	# build the package
 	export CFLAGS='-std=c++11' && dpkg-buildpackage -i -I -rfakeroot
 
-standalone: standalone.cc nmie.cc
-	export CFLAGS='-std=c++11' && c++ -DNDEBUG -O2 -Wall -std=c++11 standalone.cc nmie.cc -lm -o scattnlay
+standalone: standalone.cc nmie.cc bessel.cc
+	export CFLAGS='-std=c++11' && c++ -DNDEBUG -O2 -Wall -std=c++11 standalone.cc nmie.cc bessel.cc -lm -o scattnlay
 	mv scattnlay ../
 
 clean:

bessel.cc

@@ -34,6 +34,22 @@
 
 namespace nmie {
   namespace bessel {
+
+    void calcZeta( std::vector< std::complex<double> >& Zeta,
+		   std::vector< std::complex<double> >& dZeta,
+		   int n,  std::complex<double>z) {
+      std::vector< std::complex<double> > csj, cdj, csy, cdy;
+      int nm;
+      csphjy (n, z, nm, csj, cdj,  csy, cdy );
+      Zeta.resize(n+1);
+      dZeta.resize(n+1);
+      auto c_i = std::complex<double>(0.0,1.0);
+      for (int i = 0; i < n+1; ++i) {
+	Zeta[i]=z*(csj[i] + c_i*csy[i]);
+	dZeta[i]=z*(cdj[i] + c_i*cdy[i])+csj[i] + c_i*csy[i];
+      }
+    }  // end of calcZeta()
+
 // !*****************************************************************************80
 //  
 //  C++ port of fortran code
@@ -100,7 +116,7 @@ namespace nmie {
 	  cdy[k] = 1.0e+300;
 	}
 	csj[0] = std::complex<double>( 1.0, 0.0);
-	cdj[1] =  std::complex<double>( 0.333333333333333, 0.0);
+	cdj[1] =  std::complex<double>( 0.3333333333333333, 0.0);
 	return;
       }
       csj[0] = std::sin ( z ) / z;

bessel.h

@@ -28,6 +28,9 @@
 #include <vector>
 namespace nmie {
   namespace bessel {
+    void calcZeta( std::vector< std::complex<double> >& Zeta,
+		   std::vector< std::complex<double> >& dZeta,
+		   int n,  std::complex<double>z);
     void csphjy (int n, std::complex<double>z, int& nm,
 		 std::vector< std::complex<double> >& csj,
 		 std::vector< std::complex<double> >& cdj,

go.sh

@@ -5,17 +5,17 @@ file=compare.cc
 echo Compile with gcc
 rm -f *.bin
 rm -f ../scattnlay
-#g++ -Ofast -std=c++11 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin -static
-# g++ -Ofast -std=c++11 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay-pg.bin -static -pg
+#g++ -Ofast -std=c++11 compare.cc nmie.cc bessel.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin -static
+# g++ -Ofast -std=c++11 compare.cc nmie.cc bessel.cc nmie-wrapper.cc -lm -lrt -o scattnlay-pg.bin -static -pg
 
 #google profiler  ######## FAST!!!
 echo Uncomment next line to compile compare.cc
-# g++ -Ofast -std=c++11 $file nmie.cc nmie-old.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+# g++ -Ofast -std=c++11 $file nmie.cc bessel.cc nmie-old.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
-#  g++ -Ofast -std=c++11 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay-g.bin -ltcmalloc -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -g
+#  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 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

nmie.cc

@@ -37,6 +37,7 @@
 //                                                                                  //
 // Hereinafter all equations numbers refer to [2]                                   //
 //**********************************************************************************//
+#include "bessel.h"
 #include "nmie.h"
 #include <array>
 #include <algorithm>
@@ -560,9 +561,9 @@ namespace nmie {
     for (int i = 0; i < nmax_ + 1; ++i) {
       dPsiX[i] = D1X[i]*PsiX[i];
       dPsiMX[i] = D1MX[i]*PsiMX[i];
-      dZetaX[i] = D3X[i]*ZetaX[i];
+      //dZetaX[i] = D3X[i]*ZetaX[i];
     }
-
+    bessel::calcZeta(ZetaX, dZetaX, nmax_, x);
     an.resize(nmax_);
     bn.resize(nmax_);
     for (int i = 0; i < nmax_; i++) {
@@ -575,15 +576,15 @@ namespace nmie {
       bn[i] = Num/Denom;      
     }
     printf("dPsiX\n");
-    for (auto a: dPsiX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    for (auto a: dPsiX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
     printf("\ndPsiMX\n");
-    for (auto a: dPsiMX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    for (auto a: dPsiMX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
     printf("\nPsiX\n");
-    for (auto a: PsiX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    for (auto a: PsiX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
     printf("\nPsiMX\n");
-    for (auto a: PsiMX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    for (auto a: PsiMX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
     printf("\nZetaX\n");
-    for (auto a: ZetaX) printf("%10.5e + %10.5ei  ",std::real(a), std::imag(a));
+    for (auto a: ZetaX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
     //   It should be
     // Columns 1 through 3:
     //   3.3333e-07 - 1.0000e+03i   2.6645e-10 - 3.0000e+06i   5.1499e-07 - 1.5000e+10i
@@ -599,7 +600,7 @@ namespace nmie {
     // -2.37587e+12 + -1.35135e+26i  -3.56380e+16 + -2.02703e+30i  -6.05846e+20 + -3.44594e+34i
 
     printf("\ndZetaX\n");
-    for (auto a: dZetaX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    for (auto a: dZetaX) printf("%10.5er%+10.5ei  ",std::real(a), std::imag(a));
 
     printf("\nsize param: %g\n", x);
 

tests/bessel/test_bessel.cc

@@ -60,19 +60,15 @@ int main(int argc, char *argv[]) {
 
 
   n = 20;
-  std::complex<double> z1(1.0, 0.0);
+
+  std::complex<double> z1(0.001, 0.0);
   nmie::bessel::csphjy (n, z1, nm, csj, cdj,  csy, cdy );
   result =  csj;
 
   printf("$$$$ REAL $$$$$$ Calculate and compare against Wolfram Alpha\n");
-  printf("j(0,1) = %16.18f\n", real(result[0]));
-  printf("WA j() = 0.841470984807896506652502321630\n");
-  printf("j(1,1) = %16.18f\n", real(result[1]));
-  printf("WA j() = 0.301168678939756789251565714187322395890252640\n");
-  printf("j(1,1) = %.14g\n", real(result[10]));
-  printf("WA j() = 7.116552640047313023966191737248811458533809572 × 10^-11\n");
-  printf("j(20,1) = %.14g\n", real(result[20]));
-  printf("WA j() = 7.537795722236872993957557741584960855614358030 × 10^-26\n");
+  printf("y(0,1) = %16.18f\n", real(result[0]));
+  printf("y(1,1) = %16.18f\n", real(result[1]));
+  printf("y(2,1) = %.14g\n", real(result[2]));