#ifndef SRC_SPECIAL_FUNCTIONS_IMPL_HPP_
#define SRC_SPECIAL_FUNCTIONS_IMPL_HPP_
//******************************************************************************
//    Copyright (C) 2009-2022  Ovidio Pena <ovidio@bytesfall.com>
//    Copyright (C) 2013-2022  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 at least one of the following references:
//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
//        a multilayered sphere," Computer Physics Communications,
//        vol. 180, Nov. 2009, pp. 2348-2354.
//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
//        calculation of electromagnetic near-field for a multilayered
//        sphere," Computer Physics Communications, vol. 214, May 2017,
//        pp. 225-230.
//
//    You should have received a copy of the GNU General Public License
//    along with this program.  If not, see <http://www.gnu.org/licenses/>.
//******************************************************************************

//******************************************************************************
// This class implements the algorithm for a multilayered sphere described by:
//    [1] W. Yang, "Improved recursive algorithm for light scattering by a
//        multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp.
//        1710-1720.
//
// You can find the description of all the used equations in:
//    [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
//        a multilayered sphere," Computer Physics Communications,
//        vol. 180, Nov. 2009, pp. 2348-2354.
//    [3] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
//        calculation of electromagnetic near-field for a multilayered
//        sphere," Computer Physics Communications, vol. 214, May 2017,
//        pp. 225-230.
//
// Hereinafter all equations numbers refer to [2]
//******************************************************************************
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

#include "nmie-precision.hpp"
#include "nmie.hpp"

namespace nmie {

//******************************************************************************
// Note, that Kapteyn seems to be too optimistic (at least by 3 digits
// in some cases) for forward recurrence, see D1test with WYang_data
//******************************************************************************
template <typename FloatType>
int evalKapteynNumberOfLostSignificantDigits(const int ni,
                                             const std::complex<FloatType> zz) {
  using std::abs;
  using std::imag;
  using std::log;
  using std::real;
  using std::round;
  using std::sqrt;
  std::complex<double> z = ConvertComplex<double>(zz);
  if (abs(z) == 0.)
    return 0;
  auto n = static_cast<double>(ni);
  auto one = std::complex<double>(1, 0);
  auto lost_digits = round(  //
      (                      //
          abs(imag(z)) - log(2.) -
          n * (real(log(z / n) + sqrt(one - pow2(z / n)) -
                    log(one + sqrt(one - pow2(z / n)))  //
                    )                                   //
               )                                        //
          ) /
      log(10.)  //
  );
  return lost_digits;
}

//******************************************************************************
template <typename FloatType>
int getNStar(int nmax, std::complex<FloatType> z, const int valid_digits) {
  if (nmax == 0)
    nmax = 1;
  int nstar = nmax;
  auto z_dp = ConvertComplex<double>(z);
  if (std::abs(z_dp) == 0.)
    return nstar;
  int forwardLoss = evalKapteynNumberOfLostSignificantDigits(nmax, z_dp);
  int increment = static_cast<int>(
      std::ceil(std::max(4 * std::pow(std::abs(z_dp), 1 / 3.0), 5.0)));
  int backwardLoss = evalKapteynNumberOfLostSignificantDigits(nstar, z_dp);
  while (backwardLoss - forwardLoss < valid_digits) {
    nstar += increment;
    backwardLoss = evalKapteynNumberOfLostSignificantDigits(nstar, z_dp);
  };
  return nstar;
}

//******************************************************************************
// Custom implementation of complex cot function to avoid overflow
// if Im(z) < 0, then it evaluates cot(z) as conj(cot(conj(z)))
// see Eqs. 10-12 of [1] for details.
// [1]H. Du, Mie-Scattering Calculation, Appl. Opt. 43, 1951 (2004).
//******************************************************************************
template <typename FloatType>
std::complex<FloatType> complex_cot(const std::complex<FloatType> z) {
  auto Remx = z.real();
  auto Immx = z.imag();
  int sign =
      (Immx > 0) ? 1 : -1;  // use complex conj if needed for exp and return
  auto exp = nmm::exp(-2 * sign * Immx);
  auto tan = nmm::tan(Remx);
  auto a = tan - exp * tan;
  auto b = 1 + exp;
  auto c = -1 + exp;
  auto d = tan + exp * tan;
  auto c_one = std::complex<FloatType>(0, 1);
  return (a * c + b * d) / (pow2(c) + pow2(d)) +
         c_one * (sign * (b * c - a * d) / (pow2(c) + pow2(d)));
}

//******************************************************************************
// Forward iteration for evaluation of ratio of the Riccati–Bessel functions
//******************************************************************************
template <typename FloatType>
void evalForwardR(const std::complex<FloatType> z,
                  std::vector<std::complex<FloatType>>& r) {
  if (r.size() < 1)
    throw std::invalid_argument(
        "We need non-zero evaluations of ratio of the Riccati–Bessel "
        "functions.\n");
  // r0 = cot(z)
  //  r[0] = nmm::cos(z)/nmm::sin(z);
  r[0] = complex_cot(z);
  for (unsigned int n = 0; n < r.size() - 1; n++) {
    r[n + 1] = static_cast<FloatType>(1) /
               (                                          //
                   static_cast<FloatType>(2 * n + 1) / z  //
                   - r[n]                                 //
               );
  }
}

//******************************************************************************
// Backward iteration for evaluation of ratio of the Riccati–Bessel functions
//******************************************************************************
template <typename FloatType>
void evalBackwardR(const std::complex<FloatType> z,
                   std::vector<std::complex<FloatType>>& r) {
  if (r.size() < 1)
    throw std::invalid_argument(
        "We need non-zero evaluations of ratio of the Riccati-Bessel "
        "functions.\n");
  int nmax = r.size() - 1;
  auto tmp = static_cast<FloatType>(2 * nmax + 1) / z;
  r[nmax] = tmp;
  for (int n = nmax - 1; n >= 0; n--) {
    r[n] = -static_cast<FloatType>(1) / r[n + 1] +
           static_cast<FloatType>(2 * n + 1) / z;
  }
  r[0] = complex_cot(z);
}

//******************************************************************************
template <typename FloatType>
void convertRtoD1(const std::complex<FloatType> z,
                  std::vector<std::complex<FloatType>>& r,
                  std::vector<std::complex<FloatType>>& D1) {
  if (D1.size() > r.size())
    throw std::invalid_argument(
        "Not enough elements in array of ratio of the Riccati–Bessel functions "
        "to convert it into logarithmic derivative array.\n");
  std::complex<FloatType> D_old;
  for (unsigned int n = 0; n < D1.size(); n++) {
    D_old = D1[n];
    D1[n] = r[n] - static_cast<FloatType>(n) / z;
    //    std::cout << "D1old - D1[n] :" << D_old-D1[n] << '\n';
  }
}

//******************************************************************************
template <typename FloatType>
void evalForwardD(const std::complex<FloatType> z,
                  std::vector<std::complex<FloatType>>& D) {
  int nmax = D.size();
  FloatType one = static_cast<FloatType>(1);
  for (int n = 1; n < nmax; n++) {
    FloatType nf = static_cast<FloatType>(n);
    D[n] = one / (nf / z - D[n - 1]) - nf / z;
  }
}

//******************************************************************************
template <typename FloatType>
void evalForwardD1(const std::complex<FloatType> z,
                   std::vector<std::complex<FloatType>>& D) {
  if (D.size() < 1)
    throw std::invalid_argument("Should have a leas one element!\n");
  D[0] = nmm::cos(z) / nmm::sin(z);
  evalForwardD(z, D);
}

//  template <typename FloatType>
//  void MultiLayerMie<FloatType>::calcD1D3(const std::complex<FloatType> z,
//                               std::vector<std::complex<FloatType> > &D1,
//                               std::vector<std::complex<FloatType> > &D3) {
//
//    // if (cabs(D1[0]) > 1.0e15) {
//    //   throw std::invalid_argument("Unstable D1! Please, try to change input
//    parameters!\n");
//    // //printf("Warning: Potentially unstable D1! Please, try to change input
//    parameters!\n");
//    // }
//
//    // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
//    PsiZeta_[0] = static_cast<FloatType>(0.5)*(static_cast<FloatType>(1.0) -
//    std::complex<FloatType>(nmm::cos(2.0*z.real()), nmm::sin(2.0*z.real()))
//                 *static_cast<FloatType>(nmm::exp(-2.0*z.imag())));
//    D3[0] = std::complex<FloatType>(0.0, 1.0);
//
//    for (int n = 1; n <= nmax_; n++) {
//      PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<FloatType>(n)*z_inv - D1[n -
//      1])
//                                   *(static_cast<FloatType>(n)*z_inv - D3[n -
//                                   1]);
//      D3[n] = D1[n] + std::complex<FloatType>(0.0, 1.0)/PsiZeta_[n];
//    }
//  }
//
//
//*******************************************************************************
// This function calculates the Riccati-Bessel functions (Psi and Zeta) for a
// complex argument (z).
// Equations (20a) - (21b)
//
// Input parameters:
//   z: Complex argument to evaluate Psi and Zeta
//   nmax: Maximum number of terms to calculate Psi and Zeta
//
// Output parameters:
//   Psi, Zeta: Riccati-Bessel functions
//*******************************************************************************
//  template <typename FloatType>
//  void MultiLayerMie<FloatType>::calcPsiZeta(std::complex<FloatType> z,
//                                  std::vector<std::complex<FloatType> > &Psi,
//                                  std::vector<std::complex<FloatType> > &Zeta)
//                                  {
//
//    std::complex<FloatType> c_i(0.0, 1.0);
//    std::vector<std::complex<FloatType> > D1(nmax_ + 1), D3(nmax_ + 1);
//
//    // First, calculate the logarithmic derivatives
//    calcD1D3(z, D1, D3);
//
//    // Now, use the upward recurrence to calculate Psi and Zeta - equations
//    (20a) - (21b) Psi[0] = nmm::sin(z); Zeta[0] = nmm::sin(z) -
//    c_i*nmm::cos(z); for (int n = 1; n <= nmax_; n++) {
//      Psi[n]  =  Psi[n - 1]*(std::complex<FloatType>(n,0.0)/z - D1[n - 1]);
//      Zeta[n] = Zeta[n - 1]*(std::complex<FloatType>(n,0.0)/z - D3[n - 1]);
//    }
//  }

//**********************************************************************************
// This function calculates Pi and Tau for a given value of cos(Theta).
// Equations (26a) - (26c)
//
// Input parameters:
//   nmax_: Maximum number of terms to calculate Pi and Tau
//   nTheta: Number of scattering angles
//   Theta: Array containing all the scattering angles where the scattering
//          amplitudes will be calculated
//
// Output parameters:
//   Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) -
//   (26c)
//**********************************************************************************
//  template <typename FloatType>
//  void MultiLayerMie<FloatType>::calcPiTau(const FloatType &costheta,
//                                std::vector<FloatType> &Pi,
//                                std::vector<FloatType> &Tau) {
//
//    int i;
//    //****************************************************//
//    // Equations (26a) - (26c)                            //
//    //****************************************************//
//    // Initialize Pi and Tau
//    Pi[0] = 1.0;  // n=1
//    Tau[0] = costheta;
//    // Calculate the actual values
//    if (nmax_ > 1) {
//      Pi[1] = 3*costheta*Pi[0]; //n=2
//      Tau[1] = 2*costheta*Pi[1] - 3*Pi[0];
//      for (i = 2; i < nmax_; i++) { //n=[3..nmax_]
//        Pi[i] = ((i + i + 1)*costheta*Pi[i - 1] - (i + 1)*Pi[i - 2])/i;
//        Tau[i] = (i + 1)*costheta*Pi[i] - (i + 2)*Pi[i - 1];
//      }
//    }
//  }  // end of MultiLayerMie::calcPiTau(...)

//**********************************************************************************
// This function calculates vector spherical harmonics (eq. 4.50, p. 95 BH),
// required to calculate the near-field parameters.
//
// Input parameters:
//   Rho: Radial distance
//   Phi: Azimuthal angle
//   Theta: Polar angle
//   rn: Either the spherical Ricatti-Bessel function of first or third kind
//   Dn: Logarithmic derivative of rn
//   Pi, Tau: Angular functions Pi and Tau
//   n: Order of vector spherical harmonics
//
// Output parameters:
//   Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics
//**********************************************************************************
//  template <typename FloatType>
//  void MultiLayerMie<FloatType>::calcSpherHarm(const std::complex<FloatType>
//  Rho, const FloatType Theta, const FloatType Phi,
//                                    const std::complex<FloatType> &rn, const
//                                    std::complex<FloatType> &Dn, const
//                                    FloatType &Pi, const FloatType &Tau, const
//                                    FloatType &n,
//                                    std::vector<std::complex<FloatType> >&
//                                    Mo1n, std::vector<std::complex<FloatType>
//                                    >& Me1n,
//                                    std::vector<std::complex<FloatType> >&
//                                    No1n, std::vector<std::complex<FloatType>
//                                    >& Ne1n) {
//
//    // using eq 4.50 in BH
//    std::complex<FloatType> c_zero(0.0, 0.0);
//
//    using nmm::sin;
//    using nmm::cos;
//    Mo1n[0] = c_zero;
//    Mo1n[1] = cos(Phi)*Pi*rn/Rho;
//    Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
//
//    Me1n[0] = c_zero;
//    Me1n[1] = -sin(Phi)*Pi*rn/Rho;
//    Me1n[2] = -cos(Phi)*Tau*rn/Rho;
//
//    No1n[0] = sin(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
//    No1n[1] = sin(Phi)*Tau*Dn*rn/Rho;
//    No1n[2] = cos(Phi)*Pi*Dn*rn/Rho;
//
//    Ne1n[0] = cos(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
//    Ne1n[1] = cos(Phi)*Tau*Dn*rn/Rho;
//    Ne1n[2] = -sin(Phi)*Pi*Dn*rn/Rho;
//  }  // end of MultiLayerMie::calcSpherHarm(...)
//

//******************************************************************************
// This functions calculate the logarithmic derivatives of the Riccati-Bessel
// functions (D1 and D3) for a complex argument (z).
// Equations (16a), (16b) and (18a) - (18d)
//
// Input parameters:
//   z: Complex argument to evaluate D1 and D3
//   nmax_: Maximum number of terms to calculate D1 and D3
//
// Output parameters:
//   D1, D3: Logarithmic derivatives of the Riccati-Bessel functions
//******************************************************************************
template <typename FloatType>
void evalDownwardD1(const std::complex<FloatType> z,
                    std::vector<std::complex<FloatType>>& D1) {
  int nmax = D1.size() - 1;
  int valid_digits = 16;
#ifdef MULTI_PRECISION
  valid_digits += MULTI_PRECISION;
#endif
  int nstar = nmie::getNStar(nmax, z, valid_digits);
  D1.resize(nstar + 1);
  // Downward recurrence for D1 - equations (16a) and (16b)
  D1[nstar] = std::complex<FloatType>(0.0, 0.0);
  std::complex<FloatType> c_one(1.0, 0.0);
  const std::complex<FloatType> z_inv = std::complex<FloatType>(1.0, 0.0) / z;
  for (unsigned int n = nstar; n > 0; n--) {
    D1[n - 1] = static_cast<FloatType>(n) * z_inv -
                c_one / (D1[n] + static_cast<FloatType>(n) * z_inv);
  }
  // Use D1[0] from upward recurrence
  D1[0] = complex_cot(z);
  D1.resize(nmax + 1);
  //  printf("D1[0] = (%16.15g, %16.15g) z=(%16.15g,%16.15g)\n",
  //  D1[0].real(),D1[0].imag(),
  //         z.real(),z.imag());
}

//******************************************************************************
template <typename FloatType>
void evalUpwardD3(const std::complex<FloatType> z,
                  const std::vector<std::complex<FloatType>>& D1,
                  std::vector<std::complex<FloatType>>& D3,
                  std::vector<std::complex<FloatType>>& PsiZeta) {
  int nmax = D1.size() - 1;
  // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
  PsiZeta[0] = static_cast<FloatType>(0.5) *
               (static_cast<FloatType>(1.0) -
                std::complex<FloatType>(nmm::cos(2.0 * z.real()),
                                        nmm::sin(2.0 * z.real())) *
                    static_cast<FloatType>(nmm::exp(-2.0 * z.imag())));
  D3[0] = std::complex<FloatType>(0.0, 1.0);
  const std::complex<FloatType> z_inv = std::complex<FloatType>(1.0, 0.0) / z;
  for (int n = 1; n <= nmax; n++) {
    PsiZeta[n] = PsiZeta[n - 1] *
                 (static_cast<FloatType>(n) * z_inv - D1[n - 1]) *
                 (static_cast<FloatType>(n) * z_inv - D3[n - 1]);
    D3[n] = D1[n] + std::complex<FloatType>(0.0, 1.0) / PsiZeta[n];
  }
}

//******************************************************************************
template <typename FloatType>
std::complex<FloatType> complex_sin(const std::complex<FloatType> z) {
  auto a = z.real();
  auto b = z.imag();
  auto i = std::complex<FloatType>(0, 1);
  using nmm::cos;
  using nmm::exp;
  using nmm::sin;
  return ((cos(a) + i * sin(a)) * exp(-b) - (cos(-a) + i * sin(-a)) * exp(b)) /
         (static_cast<FloatType>(2) * i);
}

//******************************************************************************
template <typename FloatType>
void evalUpwardPsi(const std::complex<FloatType> z,
                   const std::vector<std::complex<FloatType>> D1,
                   std::vector<std::complex<FloatType>>& Psi) {
  int nmax = Psi.size() - 1;
  // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a)
  // - (21b)
  Psi[0] = complex_sin(z);
  for (int n = 1; n <= nmax; n++) {
    Psi[n] = Psi[n - 1] * (std::complex<FloatType>(n, 0.0) / z - D1[n - 1]);
  }
}

//******************************************************************************
// Sometimes in literature Zeta is also denoted as Ksi, it is a Riccati-Bessel
// function of third kind.
//******************************************************************************
template <typename FloatType>
void evalUpwardZeta(const std::complex<FloatType> z,
                    const std::vector<std::complex<FloatType>> D3,
                    std::vector<std::complex<FloatType>>& Zeta) {
  int nmax = Zeta.size() - 1;
  std::complex<FloatType> c_i(0.0, 1.0);
  // Now, use the upward recurrence to calculate Zeta and Zeta - equations (20a)
  // - (21b)
  Zeta[0] = nmm::sin(z) - c_i * nmm::cos(z);
  for (int n = 1; n <= nmax; n++) {
    Zeta[n] = Zeta[n - 1] * (std::complex<FloatType>(n, 0.0) / z - D3[n - 1]);
  }
}

//******************************************************************************
// void evalForwardRiccatiBessel(const FloatType x, const FloatType first, const
// FloatType second,
//                               std::vector<FloatType> &values) {
//   values[0] = first;
//   values[1] = second;
//   int nmax = values.size();
//   for (int i = 1; i < nmax-1; i++) {
//     values[i+1] = (1 + 2*i) * values[i]/x - values[i-1];
//   }
// }
//
//******************************************************************************
// void evalChi(const FloatType x, std::vector<FloatType> &Chi) {
//   auto first = nmm::cos(x);
//   auto second = first/x + nmm::sin(x);
//   evalForwardRiccatiBessel(x, first, second, Chi);
// }
//
//******************************************************************************
// void evalPsi(const FloatType x, std::vector<FloatType> &Psi) {
//   auto first = nmm::sin(x);
//   auto second = first/x - nmm::cos(x);
//   evalForwardRiccatiBessel(x, first, second, Psi);
// }
//
//******************************************************************************
// void composeZeta(const std::vector<FloatType> &Psi,
//                  const std::vector<FloatType> &Chi,
//                  std::vector< std::complex<FloatType>> &Zeta) {
//   int nmax = Zeta.size();
//   for (int i = 0; i < nmax; i++) {
//     Zeta[i] = std::complex<FloatType > (Psi[i], Chi[i]);
//   }
// }
//******************************************************************************
template <typename FloatType>
void evalPsiZetaD1D3(const std::complex<FloatType> cxd,
                     std::vector<std::complex<FloatType>>& Psi,
                     std::vector<std::complex<FloatType>>& Zeta,
                     std::vector<std::complex<FloatType>>& D1,
                     std::vector<std::complex<FloatType>>& D3) {
  int nmax = Psi.size() - 1;
  std::vector<std::complex<FloatType>> PsiZeta(nmax + 1);

  for (int n = 0; n <= nmax; n++) {
    D1[n] = std::complex<FloatType>(0.0, -1.0);
    D3[n] = std::complex<FloatType>(0.0, 1.0);
  }

  evalDownwardD1(cxd, D1);
  evalUpwardPsi(cxd, D1, Psi);
  evalUpwardD3(cxd, D1, D3, PsiZeta);
  for (unsigned int i = 0; i < Zeta.size(); i++) {
    Zeta[i] = PsiZeta[i] / Psi[i];
  }
}

}  // end of namespace nmie
#endif  // SRC_SPECIAL_FUNCTIONS_IMPL_HPP_