Tests for diff(Zeta)

bessel.cc

@@ -35,9 +35,8 @@
 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 calcZeta(int n,  std::complex<double>z,  std::vector< std::complex<double> >& Zeta,
+		   std::vector< std::complex<double> >& dZeta) {
       std::vector< std::complex<double> > csj, cdj, csy, cdy;
       int nm;
       csphjy (n, z, nm, csj, cdj,  csy, cdy );
@@ -50,6 +49,19 @@ namespace nmie {
       }
     }  // end of calcZeta()
 
+    void calcPsi(int n,  std::complex<double>z,  std::vector< std::complex<double> >& Psi,
+		   std::vector< std::complex<double> >& dPsi) {
+      std::vector< std::complex<double> > csj, cdj, csy, cdy;
+      int nm;
+      csphjy (n, z, nm, csj, cdj,  csy, cdy );
+      Psi.resize(n+1);
+      dPsi.resize(n+1);
+      for (int i = 0; i < n+1; ++i) {
+	Psi[i]=z*(csj[i]);
+	dPsi[i]=z*(cdj[i])+csj[i];
+      }
+    }  // end of calcPsi()
+
 // !*****************************************************************************80
 //  
 //  C++ port of fortran code

bessel.h

@@ -28,9 +28,12 @@
 #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 calcZeta(int n,  std::complex<double>z,
+		  std::vector< std::complex<double> >& Zeta,
+		  std::vector< std::complex<double> >& dZeta);
+    void calcPsi(int n,  std::complex<double>z,
+		  std::vector< std::complex<double> >& Psi,
+		  std::vector< std::complex<double> >& dPsi);
     void csphjy (int n, std::complex<double>z, int& nm,
 		 std::vector< std::complex<double> >& csj,
 		 std::vector< std::complex<double> >& cdj,

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 @@ PROGRAM='./scattnlay.bin'
 # ./$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

@@ -563,7 +563,7 @@ namespace nmie {
       dPsiMX[i] = D1MX[i]*PsiMX[i];
       //dZetaX[i] = D3X[i]*ZetaX[i];
     }
-    bessel::calcZeta(ZetaX, dZetaX, nmax_, x);
+    bessel::calcZeta(nmax_, x, ZetaX, dZetaX);
     an.resize(nmax_);
     bn.resize(nmax_);
     for (int i = 0; i < nmax_; i++) {
@@ -650,7 +650,7 @@ namespace nmie {
       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);
+    bessel::calcZeta(nmax_, z, ZetaZ, dZetaZ );
     for (int n = 0; n < nmax_+1; ++n) {
       D3[n]=dZetaZ[n]/ZetaZ[n];
     }
@@ -689,7 +689,7 @@ namespace nmie {
     }
     
     std::vector<std::complex<double> >  dZetaZ(nmax_ + 1);
-    bessel::calcZeta(Zeta, dZetaZ, nmax_, z);
+    bessel::calcZeta(nmax_, z, Zeta, dZetaZ);
 
 
   }

tests/bessel/test_bessel.cc

@@ -35,43 +35,287 @@
 
 int main(int argc, char *argv[]) {
   
-  int n = 10;
-  std::complex<double> z(1.0, 2.0);
+  int n_small = 1, n_big = 15;
+  int n = n_big+2, n_print;
+  std::complex<double> z_small(0.0012, 0.0034);
+  std::complex<double> z_medium(1.0, 2.0);
+  std::complex<double> z_big(18.0, 17.0);
+  std::complex<double> z_big_real(81.0, 0.0);
+  std::complex<double> z_big_imag(0.0, 13.0);
+  std::complex<double> z;
+
   int nm;
-  std::vector< std::complex<double> > csj, cdj, csy, cdy;
+  // csj and cdj - complex argument spherical bessel first kind and derivative
+  // csy and cdy - complex argument spherical bessel second kind and derivative
+  std::vector< std::complex<double> > csj, cdj, csy, cdy, csh;
+  std::vector<std::complex<double> > Psi, Zeta, dPsi, dZeta;
+
+  auto c_i = std::complex<double>(0.0,1.0);
+
+  printf("===== Small ========\n");
+  z = z_small;
+  printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
+  nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+  csh.resize(csj.size());
+  for (unsigned int i = 0; i < csj.size(); ++i)
+    csh[i] = csj[i]+c_i*csy[i];
+  nmie::bessel::calcPsi(n, z, Psi, dPsi);
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+
+  // n_print = n_small;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf("WA csj =         e   r        e   i;  WA csy =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf(" WA csj =         e   r        e   i;   WA csy =         e   r        e   i\n");
+
+  // n_print = n_small;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj =         e   r        e   i;\n");
+  // n_print = n_big;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj  =         e   r        e   i;\n");
+    
+  // n_print = n_small;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf("WA Psi =         e   r        e   i;  WA Zeta =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf(" WA Psi =         e   r        e   i;   WA Zeta =         e   r        e   i\n");
+
+  n_print = n_small;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
+	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf("WA dPsi =         e   r        e   i;  WA dZeta = +4.8284 e+04r-5.98822e+04i\n");
+  n_print = n_big;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
+	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf(" WA dPsi =         e   r        e   i;   WA dZeta = +8.58553e+55r+7.47397e+55i\n");
+
+
+
+
+  printf("===== Medium ========\n");
+  z = z_medium;
+  printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
+  nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+  csh.resize(csj.size());
+  for (unsigned int i = 0; i < csj.size(); ++i)
+    csh[i] = csj[i]+c_i*csy[i];
+  nmie::bessel::calcPsi(n, z, Psi, dPsi);
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+
+  // n_print = n_small;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf("WA csj =         e   r        e   i;  WA csy =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf(" WA csj =         e   r        e   i;   WA csy =         e   r        e   i\n");
+
+  // n_print = n_small;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj =         e   r        e   i;\n");
+  // n_print = n_big;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj  =         e   r        e   i;\n");
+    
+  // n_print = n_small;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf("WA Psi =         e   r        e   i;  WA Zeta =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf(" WA Psi =         e   r        e   i;   WA Zeta =         e   r        e   i\n");
+
+  n_print = n_small;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf("WA dPsi =         e   r        e   i;  WA dZeta = +1.99423e-01r-7.01483e-02i\n");
+  n_print = n_big;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf(" WA dPsi =         e   r        e   i;   WA dZeta = -2.11402e+11r+8.34528e+10i\n");
+
+
+  printf("===== Big ========\n");
+  z = z_big;
+  printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
   nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+  csh.resize(csj.size());
+  for (unsigned int i = 0; i < csj.size(); ++i)
+    csh[i] = csj[i]+c_i*csy[i];
+  nmie::bessel::calcPsi(n, z, Psi, dPsi);
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
 
-  auto result =  csj;
-  printf("===========Calculate and compare against Wolfram Alpha\n");
-  printf("j(0,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-	 real(result[0]),
-	 imag(result[0]));
-  printf("WA j() = re(1.4169961192118759)\n         im(-0.874391197002)\n");
+  // n_print = n_small;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf("WA csj =         e   r        e   i;  WA csy =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf(" WA csj =         e   r        e   i;   WA csy =         e   r        e   i\n");
 
-  printf("j(1,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-	 real(result[1]),
-	 imag(result[1]));
-  printf("WA j() = re(0.74785726329830368)\n         im(0.68179207555304)\n");
+  // n_print = n_small;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj =         e   r        e   i;\n");
+  // n_print = n_big;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj  =         e   r        e   i;\n");
+    
+  // n_print = n_small;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf("WA Psi =         e   r        e   i;  WA Zeta =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf(" WA Psi =         e   r        e   i;   WA Zeta =         e   r        e   i\n");
 
-  printf("j(4,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-	 real(result[4]),
-	 imag(result[4]));
-  printf("WA j() = re(-0.01352410550046)\n         im(-0.027169663050653)\n");
+  n_print = n_small;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf("WA dPsi =         e   r        e   i;  WA dZeta = -3.11025e-08r-2.90558e-08i\n");
+  n_print = n_big;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf(" WA dPsi =         e   r        e   i;   WA dZeta = -6.13497e-07r-1.14967e-06i\n");
 
 
-  n = 20;
+  printf("===== Big real========\n");
+  z = z_big_real;
+  printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
+  nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+  csh.resize(csj.size());
+  for (unsigned int i = 0; i < csj.size(); ++i)
+    csh[i] = csj[i]+c_i*csy[i];
+  nmie::bessel::calcPsi(n, z, Psi, dPsi);
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
+
+  // n_print = n_small;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf("WA csj =         e   r        e   i;  WA csy =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf(" WA csj =         e   r        e   i;   WA csy =         e   r        e   i\n");
+
+  // n_print = n_small;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj =         e   r        e   i;\n");
+  // n_print = n_big;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj  =         e   r        e   i;\n");
+    
+  // n_print = n_small;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf("WA Psi =         e   r        e   i;  WA Zeta =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf(" WA Psi =         e   r        e   i;   WA Zeta =         e   r        e   i\n");
 
-  std::complex<double> z1(0.001, 0.0);
-  nmie::bessel::csphjy (n, z1, nm, csj, cdj,  csy, cdy );
-  result =  csj;
+  n_print = n_small;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf("WA dPsi =         e   r        e   i;  WA dZeta = -6.20203e-01r-7.84344e-01i\n");
+  n_print = n_big;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf(" WA dPsi =         e   r        e   i;   WA dZeta = -7.13982e-01r+6.86858e-01i\n");
 
-  printf("$$$$ REAL $$$$$$ Calculate and compare against Wolfram Alpha\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]));
+  printf("===== Big imag========\n");
+  z = z_big_imag;
+  printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
+  nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+  csh.resize(csj.size());
+  for (unsigned int i = 0; i < csj.size(); ++i)
+    csh[i] = csj[i]+c_i*csy[i];
+  nmie::bessel::calcPsi(n, z, Psi, dPsi);
+  nmie::bessel::calcZeta(n, z, Zeta, dZeta);  
 
+  // n_print = n_small;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf("WA csj =         e   r        e   i;  WA csy =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("csj[%i] = %+10.5er%+10.5ei;  csy[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 csj[n_print].real(), csj[n_print].imag(), n_print,
+  // 	 csy[n_print].real(), csy[n_print].imag());
+  // printf(" WA csj =         e   r        e   i;   WA csy =         e   r        e   i\n");
 
+  // n_print = n_small;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj =         e   r        e   i;\n");
+  // n_print = n_big;
+  // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 csh[n_print].real(), csh[n_print].imag());
+  // printf("WA csj  =         e   r        e   i;\n");
     
+  // n_print = n_small;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf("WA Psi =         e   r        e   i;  WA Zeta =         e   r        e   i\n");
+  // n_print = n_big;
+  // printf("Psi[%i] = %+10.5er%+10.5ei;  Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  // 	 Psi[n_print].real(), Psi[n_print].imag(), n_print,
+  // 	 Zeta[n_print].real(), Zeta[n_print].imag());
+  // printf(" WA Psi =         e   r        e   i;   WA Zeta =         e   r        e   i\n");
+
+  n_print = n_small;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf("WA dPsi =         e   r        e   i;  WA dZeta = -2.47233e-22r-2.44758e-06i\n");
+  n_print = n_big;
+  printf("dPsi[%i] = %+10.5er%+10.5ei;  dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
+  	 dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
+  	 dZeta[n_print].real(), dZeta[n_print].imag());
+  printf(" WA dPsi =         e   r        e   i;   WA dZeta =  8.39449e-19r+1.28767e-02i\n");
 
 }