/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

// Third party includes
#include <boost/math/constants/constants.hpp>

// VOTCA includes
#include <votca/tools/constants.h>

// Local VOTCA includes
#include <votca/xtp/orbitals.h>

// Local private VOTCA includes
#include "spectrum.h"

namespace votca {
namespace xtp {

void Spectrum::ParseOptions(const tools::Property& options) {

  // orbitals file or pure DFT output
  _orbfile = options.ifExistsReturnElseReturnDefault<std::string>(
      ".orbitals", _job_name + ".orb");

  _output_file = options.ifExistsReturnElseReturnDefault<std::string>(
      ".output", _job_name + "_spectrum.dat");

  _n_pt = options.get(".points").as<Index>();
  _lower = options.get(".lower").as<double>();
  _upper = options.get(".upper").as<double>();
  _fwhm = options.get(".fwhm").as<double>();

  _spectrum_type = options.get(".type").as<std::string>();
  _minexc = options.get(".minexc").as<Index>();
  _maxexc = options.get(".maxexc").as<Index>();
  _shiftby = options.get(".shift").as<double>();
}

bool Spectrum::Run() {
  _log.setReportLevel(Log::current_level);
  _log.setMultithreading(true);

  _log.setCommonPreface("\n... ...");

  XTP_LOG(Log::error, _log)
      << "Calculating absorption spectrum plot " << _orbfile << std::flush;

  Orbitals orbitals;
  // load the QM data from serialized orbitals object
  XTP_LOG(Log::error, _log)
      << " Loading QM data from " << _orbfile << std::flush;
  orbitals.ReadFromCpt(_orbfile);

  // check if orbitals contains singlet energies and transition dipoles
  if (!orbitals.hasBSESinglets()) {
    throw std::runtime_error(
        "BSE singlet energies not stored in QM data file!");
  }

  if (!orbitals.hasTransitionDipoles()) {
    throw std::runtime_error(
        "BSE transition dipoles not stored in QM data file!");
  }

  const Eigen::VectorXd BSESingletEnergies =
      orbitals.BSESinglets().eigenvalues() * tools::conv::hrt2ev;
  const std::vector<Eigen::Vector3d>& TransitionDipoles =
      orbitals.TransitionDipoles();
  Eigen::VectorXd osc = orbitals.Oscillatorstrengths();

  if (_maxexc > Index(TransitionDipoles.size())) {
    _maxexc = Index(TransitionDipoles.size()) - 1;
  }

  Index n_exc = _maxexc - _minexc + 1;
  XTP_LOG(Log::error, _log)
      << " Considering " << n_exc << " excitation with max energy "
      << BSESingletEnergies(_maxexc) << " eV / min wave length "
      << evtonm(BSESingletEnergies[_maxexc - 1]) << " nm" << std::flush;

  /*
   *
   * For a single excitation, broaden by Lineshape function L(v-W)
   *    eps(v) = f * L(v-W)
   *
   * where
   *       v: energy
   *       f: oscillator strength in dipole-length gauge
   *       W: excitation energy
   *
   * Lineshape function depend on FWHM and can be
   *
   *      Gaussian
   *          L(v-W) = 1/(sqrt(2pi)sigma) * exp(-0.5 (v-W)^2/sigma^2
   *
   *
   *            with sigma: derived from FWHM (FWHM/2.3548)
   *
   *     Lorentzian
   *          L(v-W) = 1/pi * 0.5 FWHM/( (v-w)^2 + 0.25*FWHM^2 )
   *
   * Full spectrum is superposition of individual spectra.
   *
   *  Alternatively, one can calculate the imaginary part of the
   *  frequency-dependent dielectric function
   *
   *   IM(eps(v)) ~ 1/v^2 * W^2 * |td|^2 * L(v-W)
   *              = 1/v^2 * W   * f      * L(v-W)
   *
   *
   */

  std::ofstream ofs(_output_file, std::ofstream::out);

  if (_spectrum_type == "energy") {
    ofs << "# E(eV)    epsGaussian    IM(eps)Gaussian   epsLorentz    "
           "Im(esp)Lorentz\n";
    for (Index i_pt = 0; i_pt <= _n_pt; i_pt++) {

      double e = (_lower + double(i_pt) * (_upper - _lower) / double(_n_pt));

      double eps_Gaussian = 0.0;
      double imeps_Gaussian = 0.0;
      double eps_Lorentzian = 0.0;
      double imeps_Lorentzian = 0.0;

      for (Index i_exc = _minexc; i_exc <= _maxexc; i_exc++) {
        eps_Gaussian +=
            osc[i_exc] *
            Gaussian(e, BSESingletEnergies(i_exc) + _shiftby, _fwhm);
        imeps_Gaussian += osc[i_exc] * BSESingletEnergies(i_exc) *
                          Gaussian(e, BSESingletEnergies(i_exc), _fwhm);
        eps_Lorentzian +=
            osc[i_exc] * Lorentzian(e, BSESingletEnergies(i_exc), _fwhm);
        imeps_Lorentzian += osc[i_exc] * BSESingletEnergies(i_exc) *
                            Lorentzian(e, BSESingletEnergies(i_exc), _fwhm);
      }

      ofs << e << "    " << eps_Gaussian << "   " << imeps_Gaussian << "   "
          << eps_Lorentzian << "   " << imeps_Lorentzian << std::endl;
    }

    XTP_LOG(Log::error, _log)
        << " Spectrum in energy range from  " << _lower << " to " << _upper
        << " eV and with broadening of FWHM " << _fwhm
        << " eV written to file  " << _output_file << std::flush;
  }

  if (_spectrum_type == "wavelength") {

    ofs << "# lambda(nm)    epsGaussian    IM(eps)Gaussian   epsLorentz    "
           "Im(esp)Lorentz\n";
    for (Index i_pt = 0; i_pt <= _n_pt; i_pt++) {

      double lambda =
          (_lower + double(i_pt) * (_upper - _lower) / double(_n_pt));
      double eps_Gaussian = 0.0;
      double imeps_Gaussian = 0.0;
      double eps_Lorentzian = 0.0;
      double imeps_Lorentzian = 0.0;

      for (Index i_exc = _minexc; i_exc <= _maxexc; i_exc++) {
        double exc_lambda = nmtoev(BSESingletEnergies(i_exc) + _shiftby);
        eps_Gaussian += osc[i_exc] * Gaussian(lambda, exc_lambda, _fwhm);
        imeps_Gaussian +=
            osc[i_exc] * exc_lambda * Gaussian(lambda, exc_lambda, _fwhm);
        eps_Lorentzian += osc[i_exc] * Lorentzian(lambda, exc_lambda, _fwhm);
        imeps_Lorentzian +=
            osc[i_exc] * exc_lambda * Lorentzian(lambda, exc_lambda, _fwhm);
      }

      ofs << lambda << "    " << eps_Gaussian << "   " << imeps_Gaussian
          << "   " << eps_Lorentzian << "   " << imeps_Lorentzian << std::endl;
    }
    XTP_LOG(Log::error, _log)
        << " Spectrum in wavelength range from  " << _lower << " to " << _upper
        << " nm and with broadening of FWHM " << _fwhm
        << " nm written to file  " << _output_file << std::flush;
  }

  ofs.close();
  return true;
}

double Spectrum::Lorentzian(double x, double center, double fwhm) {
  return 0.5 * fwhm / (std::pow(x - center, 2) + 0.25 * fwhm * fwhm) /
         boost::math::constants::pi<double>();
}

double Spectrum::Gaussian(double x, double center, double fwhm) {
  // FWHM = 2*sqrt(2 ln2) sigma = 2.3548 sigma
  double sigma = fwhm / 2.3548;
  return std::exp(-0.5 * std::pow((x - center) / sigma, 2)) / sigma /
         sqrt(2.0 * boost::math::constants::pi<double>());
}

double Spectrum::evtonm(double eV) { return 1241.0 / eV; }

double Spectrum::evtoinvcm(double eV) { return 8065.73 * eV; }

double Spectrum::nmtoinvcm(double nm) { return 1241.0 * 8065.73 / nm; }

double Spectrum::invcmtonm(double invcm) { return 1.0e7 / invcm; }

double Spectrum::nmtoev(double nm) { return 1241.0 / nm; }

}  // namespace xtp
}  // namespace votca