/*! Type II quasisolitons in Lieb-Liniger: velocity_lifetime

  Purpose:
  Compute the expected group velocity, and lifetime from time series
  and produces the .vgtau file containing this information
  for (multi-)type II hole wavepackets in Lieb-Liniger.

  This executable requires Abacus version 2.

  See README for compilation instructions.

  Copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.

  This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero 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 Affero General Public License for more details.

  You should have received a copy of the GNU Affero General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
*/


import std;

import abacus;


int main(int argc, char* argv[])
{
  using namespace std::complex_literals;

  if (argc != 8) {
    std::cout << "Executable velocity_lifetime\n"
	      << " part of the Type II quasisolitons in Lieb-Liniger software suite\n"
	      << " copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.\n";
    std::cout << "\nPurpose:\n"
	      << " Produce the .vgtau file containing group velocity,\n"
	      << " time series coefficients and resulting lifetime estimates\n"
	      << " for (multi-)type II hole wavepackets in Lieb-Liniger.\n";

    std::cout << "\nPrerequirements:\n----------------\n"
	      << " - .states and .rho files produced by executable rho\n"
	      << " - .amplitudes file for the required protocol, produced by executable amplitudes_[protocol type]\n";

    std::cout << "\nUsage:\n------\n";
    std::cout << "velocity_lifetime <c> <L> <N> <nr holes> <width> <offset> <protocol>\n\n";

    int warg { 16 }, wtype { 10 }, wcons { 32 };
    std::cout << std::left << std::setw(warg) << "Argument" << std::setw(wtype) << "Type"
	      << std::setw(wcons) << "Constraints" << "Description\n";
    std::cout << std::left << std::setw(warg) << "--------" << std::setw(wtype) << "----"
	      << std::setw(wcons) << "-----------" << "-----------\n";
    std::cout << std::left << std::setw(warg) << "c" << std::setw(wtype) << "Real"
	      << std::setw(wcons) << "> 0"
	      << "Value of the interaction parameter\n";
    std::cout << std::left << std::setw(warg) << "L" << std::setw(wtype) << "Real"
	      << std::setw(wcons) << "> 0"
	      << "System size\n";
    std::cout << std::left << std::setw(warg) << "N" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "> 0"
	      << "Number of particles\n";
    std::cout << std::left << std::setw(warg) << "nr holes" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "1 <= nr holes <= N"
	      << "Number of holes (Type II modes)\n";
    std::cout << std::left << std::setw(warg) << "width" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "nr holes < width <= N"
	      << "Width of the hole window\n";
    std::cout << std::left << std::setw(warg) << "offset" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "0 <= offset <= N-width"
	      << "Offset of the hole window w/r to the right Fermi edge\n";
    std::cout << std::left << std::setw(warg) << "protocol" << std::setw(wtype) << "string"
	      << std::setw(wcons) << "(see amplitudes executables)"
	      << "Protocol used for defining the amplitudes in the wavepacket\n";

    return 0;
  }

  std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  Real c_ { std::stold(argv[1]) };
  Real L_ { std::stold(argv[2]) };
  int N_ { std::stoi(argv[3]) };
  int nholes_ { std::stoi(argv[4]) };
  int width_ { std::stoi(argv[5]) };
  int offset_ { std::stoi(argv[6]) };
  std::string protocol_ { argv[7] };   // protocol used to define the amplitudes


  int nr_states_;
  std::vector<std::string> label_;
  std::vector<int> iK_;
  std::vector<Real> E_;

  std::vector<std::complex<Real>> amplitude;

  std::map<std::string, std::map<std::string, std::complex<Real>>> rho_ME_;


  std::stringstream filename_base;
  filename_base << "c_" << c_ << "_N_" << N_ << "_L_" << L_ << "_nholes_" << nholes_;
  filename_base << "_width_" << width_;
  filename_base << "_offset_" << offset_;

  std::stringstream states_filename;
  states_filename << filename_base.str() << ".states";
  std::ifstream states_file;
  states_file.open(states_filename.str());
  states_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  std::stringstream amplitudes_filename;
  amplitudes_filename << filename_base.str() << "_" << protocol_ << ".amplitudes";
  std::ifstream amplitudes_file;
  amplitudes_file.open(amplitudes_filename.str());
  amplitudes_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  std::stringstream rho_ME_filename;
  rho_ME_filename << filename_base.str() << ".rho";
  std::ifstream rho_ME_file;
  rho_ME_file.open(rho_ME_filename.str());
  rho_ME_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);


  // Input the states info
  std::string tmp_label;
  int tmp_iK;
  Real tmp_E;

  states_file >> tmp_label;
  do {
    states_file >> tmp_iK >> tmp_E;
    label_.push_back(tmp_label);
    iK_.push_back(tmp_iK);
    E_.push_back(tmp_E);
  } while (states_file >> tmp_label);
  states_file.close();

  nr_states_ = int(label_.size());


  // Input the amplitudes

  std::complex<Real> tmp_amplitude;
  for (int is { 0 }; is < nr_states_; ++is) {
    amplitudes_file >> tmp_amplitude;
    amplitude.push_back(tmp_amplitude);
  }
  amplitudes_file.close();


  // Input the density operator matrix elements

  for (int ibra { 0 }; ibra < nr_states_; ++ibra) {
    for (int iket { 0 }; iket <= ibra; ++iket) {
      rho_ME_file >> rho_ME_[label_[ibra]][label_[iket]];
    }
  }
  rho_ME_file.close();


  // Compute the weighed averages of energy*k and momentum^2
  // Logic: in the time power expansion, we put
  // (d/dt)(d/dx) rho(x+vgt, t)
  // to zero.

  std::complex<Real> omegakBar { 0 };
  std::complex<Real> ksqBar { 0 };

  for (int ibra { 0 }; ibra < nr_states_; ++ibra) {
    for (int iket { 0 }; iket < ibra; ++iket) {

      omegakBar += std::conj(amplitude[ibra]) * amplitude[iket]
	* rho_ME_[label_[ibra]][label_[iket]]
	* (E_[ibra] - E_[iket]) * (twopi_r * (iK_[ibra] - iK_[iket])/L_);

      ksqBar += std::conj(amplitude[ibra]) * amplitude[iket]
	* rho_ME_[label_[ibra]][label_[iket]]
	* std::pow(twopi_r * (iK_[ibra] - iK_[iket])/L_, 2);
    }
  }


  // Define the group velocity as vg putting linear in t term to zero

  Real vg { std::real(omegakBar)/std::real(ksqBar) };

  // Compute the weighed average of moments of Galilean-shifted energy

  int maxpower { 16 };
  Real omegaminvgk;
  std::vector<std::complex<Real>> moment(maxpower+1, std::complex<Real>(0));

  for (int ibra { 0 }; ibra < nr_states_; ++ibra) {
    for (int iket { 0 }; iket < ibra; ++iket) {

      omegaminvgk = E_[ibra] - E_[iket] - vg * twopi_r * (iK_[ibra] - iK_[iket])/L_;

      for (int power { 0 }; power <= maxpower; power++) {
	moment[power] += std::conj(amplitude[ibra]) * amplitude[iket]
	  * rho_ME_[label_[ibra]][label_[iket]]
	  * std::pow(omegaminvgk, power);
      }
    }
  }


  // Extract rescaled coefficients of the power series in t

  std::vector<Real> coefficients(maxpower+1);
  for (int power { 0 }; power <= maxpower; power++) {
    coefficients[power] = 2* std::real(std::pow(1_ir, power) * moment[power])/std::tgamma(power+1);
  }

  // Estimated radius of convergence:
  // time at which order 2, 4, 6 terms are ~ 0.1 (forgetting about higher terms)
  Real t10c2 { std::pow(std::abs(0.1/coefficients[2]), 0.5) };
  Real t10c4 { std::pow(std::abs(0.1/coefficients[4]), 0.25) };
  Real t10c6 { std::pow(std::abs(0.1/coefficients[6]), 1.0/6) };
  // time at which order 2, 4, 6 terms are ~ 0.1 of original depletion depth, which is |coefficients[0]|
  Real t10c2r { std::pow(std::abs(0.1*coefficients[0]/coefficients[2]), 0.5) };
  Real t10c4r { std::pow(std::abs(0.1*coefficients[0]/coefficients[4]), 0.25) };
  Real t10c6r { std::pow(std::abs(0.1*coefficients[0]/coefficients[6]), 1.0/6) };

  std::cout << std::setprecision(5);
  std::cout << "amplitudes+matrix elements-aware group velocity  vg = " << vg << std::endl;
  std::cout << "quadratic time dependence coefficient c2 = " << coefficients[2] << std::endl;
  std::cout << "times at which a given order term is 0.1: t=t10c2 when c2 * t^2 = 0.1 (and similarly for t10c4, t10c6)" << std::endl;
  std::cout << "t10c2 (t at which o2 ~ 0.1) = " << t10c2 << std::endl;
  std::cout << "t10c4 (t at which o4 ~ 0.1) = " << t10c4 << std::endl;
  std::cout << "t10c6 (t at which o6 ~ 0.1) = " << t10c6 << std::endl;
  std::cout << "t10c2r (t at which o2 ~ 0.1*o0) = " << t10c2r << std::endl;
  std::cout << "t10c4r (t at which o4 ~ 0.1*o0) = " << t10c4r << std::endl;
  std::cout << "t10c6r (t at which o6 ~ 0.1*o0) = " << t10c6r << std::endl;
  std::cout << coefficients << std::endl;


  // Define the output file

  std::stringstream vgtau_filename;
  vgtau_filename << filename_base.str() << "_" << protocol_ << ".vgtau";
  std::ofstream vgtau_file;
  vgtau_file.open(vgtau_filename.str(), std::ios::out | std::ios::trunc);
  vgtau_file << std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  vgtau_file << "vg\t" << vg << std::endl;
  vgtau_file << "t10c2\t" << t10c2 << std::endl;
  vgtau_file << "t10c4\t" << t10c4 << std::endl;
  vgtau_file << "t10c6\t" << t10c6 << std::endl;
  vgtau_file << "t10c2r\t" << t10c2r << std::endl;
  vgtau_file << "t10c4r\t" << t10c4r << std::endl;
  vgtau_file << "t10c6r\t" << t10c6r << std::endl;
  for (int power { 0 }; power <= maxpower; power++) {
    vgtau_file << std::endl << "c" << power << "\t" << std::real(coefficients[power]);
  }

  vgtau_file.close();

  return 0;
}