//**********************************************************************************//
//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
//                                                                                  //
//    This file is part of scattnlay                                                //
//                                                                                  //
//    This program is free software: you can redistribute it and/or modify          //
//    it under the terms of the GNU General Public License as published by          //
//    the Free Software Foundation, either version 3 of the License, or             //
//    (at your option) any later version.                                           //
//                                                                                  //
//    This program is distributed in the hope that it will be useful,               //
//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
//    GNU General Public License for more details.                                  //
//                                                                                  //
//    The only additional remark is that we expect that all publications            //
//    describing work using this software, or all commercial products               //
//    using it, cite the following reference:                                       //
//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
//        a multilayered sphere," Computer Physics Communications,                  //
//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
//                                                                                  //
//    You should have received a copy of the GNU General Public License             //
//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
//**********************************************************************************//
#include <complex>
#include <cmath>
#include <cstdio>
#include <stdexcept>
#include <vector>

#include "../../bessel.h"
#include "../../nmie.h"

int nmax_ = -1;
//Temporary variables
std::vector<std::complex<double> > PsiZeta_;

void calcD1D3(const std::complex<double> z, std::vector<std::complex<double> > &D1, std::vector<std::complex<double> > &D3) {
    D1[nmax_] = std::complex<double>(0.0, 0.0);
    const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
    for (int n = nmax_; n > 0; n--) {
      D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
    }
    if (std::abs(D1[0]) > 100000.0)
      throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");

    PsiZeta_.resize(nmax_ + 1);

    // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
    PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
                      *std::exp(-2.0*z.imag()));
    D3[0] = std::complex<double>(0.0, 1.0);
    for (int n = 1; n <= nmax_; n++) {
      PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
                                   *(static_cast<double>(n)*zinv - D3[n - 1]);
      D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
    }
  }

  void calcPsiZeta(std::complex<double> z, std::vector<std::complex<double> > &Psi, std::vector<std::complex<double> > &Zeta) {
    std::complex<double> c_i(0.0, 1.0);
    std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
    calcD1D3(z, D1, D3);
    Psi.resize(nmax_+1);
    Zeta.resize(nmax_+1);

    Psi[0] = std::sin(z);
    Zeta[0] = std::sin(z) - c_i*std::cos(z);
    for (int n = 1; n <= nmax_; n++) {
      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]);
    }
  }


int main(int argc, char *argv[]) {

  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;

  std::vector<std::complex<double> > D1, D3;

  int nm;
  // 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 = 0.800005e-03r+0.226667e-02i;  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 = +1.75388e-53r-6.94429e-54i;   WA dZeta = +8.58553e+55r+7.47397e+55i\n");

  nmax_ = nm;
  printf("----- Scattnlay (nmax_=%d)-----\n",nmax_);
  calcPsiZeta(z, Psi, Zeta);

  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("===== 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 = 2.41792e   r+1.27781 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 = -1.03313e-11r-1.13267e-11i;   WA dZeta = -2.11402e+11r+8.34528e+10i\n");

  nmax_ = nm;
  printf("----- Scattnlay (nmax_=%d)-----\n",nmax_);
  calcPsiZeta(z, Psi, Zeta);

  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("===== 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);

  // 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 = -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 = -2.47740e+05r+3.05899e+05i;   WA dZeta = -6.13497e-07r-1.14967e-06i\n");

  nmax_ = 100;
  printf("----- Scattnlay (nmax_=%d)-----\n",nmax_);
  calcPsiZeta(z, Psi, Zeta);

  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("===== 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());
  //
  //SphericalHankelH1(1,81)*81
  printf("WA Psi =-0.784462377 e        e   i;  WA Zeta = -0.784462377r+0.620299278i\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 = 0.699946236744       e   i;   WA Zeta = 0.69994623  r+0.727239996i\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 = -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");

  nmax_ = 100 ;
  printf("----- Scattnlay (nmax_=%d)-----\n",nmax_);
  calcPsiZeta(z, Psi, Zeta);

  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("===== 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");

}