/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  ABACUS_LiebLin.h

Purpose:  Declares LiebLin gas-related classes and functions.

***********************************************************/

#ifndef ABACUS_LIEBLIN_H
#define ABACUS_LIEBLIN_H

#include "ABACUS.h"

namespace ABACUS {

  // First, some global constants...

  const DP ITER_REQ_PREC_LIEBLIN = 1.0e+4 * MACHINE_EPS_SQ;

  const int LIEBLIN_Ix2_MIN = -1000000; // Like a UV cutoff. Assumption: never reached in scanning.
  const int LIEBLIN_Ix2_MAX = -LIEBLIN_Ix2_MIN;


  //***********************************************************************

  /**
     Eigenstates of the Lieb-Liniger model in the repulsive \f$(c > 0)\f$ regime.

     For numerical convenience, rapidities are rescaled by the interaction parameter \f$c\f$.
     In equations throughout this documentation, \f$\tilde{\lambda} \equiv \lambda/c\f$.
     In the code, these rescaled rapidities are denoted `lambdaoc` (read: "lambda over c").
   */

  class LiebLin_Bethe_State {

  public:
    DP c_int;  // interaction parameter
    DP L;
    DP cxL;
    int N;
    std::string label;
    Vect<int> Ix2_available; // quantum numbers which are allowable but not occupied
    Vect<int> index_first_hole_to_right;
    Vect<int> displacement;
    Vect<int> Ix2;
    Vect<DP> lambdaoc;
    Vect<DP> S; // scattering sum
    Vect<DP> dSdlambdaoc; // its derivative
    DP diffsq;
    DP prec;
    int conv;
    int iter_Newton;
    DP E;
    int iK;
    DP K;
    DP lnnorm;

  public:
    LiebLin_Bethe_State ();
    LiebLin_Bethe_State (DP c_int_ref, DP L_ref, int N_ref);
    LiebLin_Bethe_State& operator= (const LiebLin_Bethe_State& RefState);

  public:
    int Charge () { return(N); };
    void Set_to_Label (std::string label_ref, const Vect<int>& OriginStateIx2);
    void Set_to_Label (std::string label_ref); // assumes OriginState == GroundState
    void Set_Label_from_Ix2 (const Vect<int>& OriginStateIx2);
    void Set_Label_Internals_from_Ix2 (const Vect<int>& OriginStateIx2);
    bool Check_Admissibility(char whichDSF); // always returns true
    void Find_Rapidities (bool reset_rapidities);
    bool Check_Rapidities();  // checks that all rapidities are not nan
    DP String_delta(); // trivially returns 0; exists to mirror spin chain function
    bool Check_Symmetry ();  // checks whether set of quantum numbers obeys { I } = { -I }
    void Compute_lnnorm ();
    void Compute_All (bool reset_rapidities);     // solves BAE, computes E, K and lnnorm
    void Set_Free_lambdaocs();
    void Iterate_BAE(DP damping);
    void Iterate_BAE_S(DP damping);
    void Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx);
    void Compute_Energy ();
    void Compute_Momentum ();
    DP Kernel (int a, int b);
    DP Kernel (DP lambdaoc_ref);
    void Build_Reduced_Gaudin_Matrix (SQMat<DP>& Gaudin_Red);
    void Build_Reduced_BEC_Quench_Gaudin_Matrix (SQMat<DP>& Gaudin_Red);
    void Annihilate_ph_pair (int ipart, int ihole, const Vect<int>& OriginStateIx2);
    void Parity_Flip ();  // takes all lambdaoc to -lambdaoc
    inline bool Set_to_Inner_Skeleton (int iKneeded, const Vect<int>& OriginIx2)
    {
      if (N < 3) ABACUSerror("N<3 incompatible with fixed momentum scanning");
      Ix2[0] = Ix2[1] - 2;
      Ix2[N-1] = Ix2[N-2] + 2;
      (*this).Compute_Momentum();
      if (iKneeded >= iK) Ix2[N-1] += 2*(iKneeded - iK);
      else Ix2[0] += 2*(iKneeded - iK);
      if (Ix2[0] < LIEBLIN_Ix2_MIN || Ix2[N-1] > LIEBLIN_Ix2_MAX) return(false);
      (*this).Set_Label_from_Ix2 (OriginIx2);
      return(true);
    }
    void Set_to_Outer_Skeleton (const Vect<int>& OriginIx2)
    {
      Ix2[0] = LIEBLIN_Ix2_MIN + (N % 2) + 1;
      Ix2[N-1] = LIEBLIN_Ix2_MAX - (N % 2) - 1;
      (*this).Set_Label_from_Ix2 (OriginIx2);
    };
  };

  inline bool Is_Inner_Skeleton (LiebLin_Bethe_State& State) {
    return (State.N >= 2 && (State.Ix2[0] == State.Ix2[1] - 2
			     || State.Ix2[State.N-1] == State.Ix2[State.N-2] + 2));
  };
  inline bool Is_Outer_Skeleton (LiebLin_Bethe_State& State) {
    return (State.N >= 2 && State.Ix2[0] == LIEBLIN_Ix2_MIN + (State.N % 2) + 1
	    && State.Ix2[State.N-1] == LIEBLIN_Ix2_MAX - (State.N % 2) - 1);
  };
  inline bool Is_Outer_Skeleton (const LiebLin_Bethe_State& State) {
    return (State.N >= 2 && State.Ix2[0] == LIEBLIN_Ix2_MIN + (State.N % 2) + 1
	    && State.Ix2[State.N-1] == LIEBLIN_Ix2_MAX - (State.N % 2) - 1);
  };

  inline bool Force_Descent (char whichDSF, LiebLin_Bethe_State& ScanState,
			     LiebLin_Bethe_State& RefState,
			     int desc_type_required, int iKmod, DP Chem_Pot)
  {
    bool forcedesc = false;

    // Force descent if we're computing density-density, we're at zero momentum
    // and we're descending with momentum preserved:
    if (whichDSF == 'd' && ScanState.iK == RefState.iK && desc_type_required > 8) forcedesc = true;

    // For BEC to c > 0 quench, g2(x=0): force first step
    else if (whichDSF == 'B' && ScanState.label == RefState.label) forcedesc = true;
    else if (whichDSF == 'C' && ScanState.label == RefState.label) forcedesc = true;

    return(forcedesc);
  }

  std::ostream& operator<< (std::ostream& s, const LiebLin_Bethe_State& state);

  // FUNCTION DECLARATIONS:

  DP Chemical_Potential (LiebLin_Bethe_State& RefState);
  DP Sumrule_Factor (char whichDSF, LiebLin_Bethe_State& RefState, DP Chem_Pot,
		     int iKmin, int iKmax);
  void Evaluate_F_Sumrule (char whichDSF, const LiebLin_Bethe_State& RefState, DP Chem_Pot,
			   int iKmin, int iKmax, const char* RAW_Cstr, const char* FSR_Cstr);
  void Evaluate_F_Sumrule (std::string prefix, char whichDSF, const LiebLin_Bethe_State& RefState,
			   DP Chem_Pot, int iKmin, int iKmax);
  void Evaluate_F_Sumrule (char whichDSF, DP c_int, DP L, int N, DP kBT, int nstates_req,
			   DP Chem_Pot, int iKmin, int iKmax, const char* FSR_Cstr);

  // in LiebLin_Utils.cc
  DP LiebLin_dE0_dc (DP c_int, DP L, int N);
  DP LiebLin_vs (DP c_int, DP L, int N);
  DP LiebLin_Dressed_Charge_N (DP c_int, DP L, int N);
  int Momentum_Right_Excitations (LiebLin_Bethe_State& ScanState);
  int Momentum_Left_Excitations (LiebLin_Bethe_State& ScanState);
  DP ln_Overlap_with_BEC (LiebLin_Bethe_State& lambda);

  DP Particle_Hole_Excitation_Cost (char whichDSF, LiebLin_Bethe_State& AveragingState);

  std::complex<DP> ln_Density_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate);
  std::complex<DP> ln_Psi_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate);
  std::complex<DP> ln_g2_ME (LiebLin_Bethe_State& mu, LiebLin_Bethe_State& lambda);

  DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax,
				     LiebLin_Bethe_State& LeftState,
				     LiebLin_Bethe_State& RefState, DP Chem_Pot,
				     std::stringstream& DAT_outfile);

  DP LiebLin_Twisted_lnnorm (Vect<std::complex<DP> >& lambdaoc, double cxL);
  std::complex<DP> LiebLin_Twisted_ln_Overlap (DP expbeta, Vect<DP> lstate_lambdaoc,
					       DP lstate_lnnorm,
					       LiebLin_Bethe_State& rstate);
  std::complex<DP> LiebLin_Twisted_ln_Overlap (std::complex<DP> expbeta,
					       Vect<std::complex<DP> > lstate_lambdaoc,
					       DP lstate_lnnorm,
					       LiebLin_Bethe_State& rstate);
  std::complex<DP> LiebLin_ln_Overlap (Vect<DP> lstate_lambdaoc, DP lstate_lnnorm,
				       LiebLin_Bethe_State& rstate);
  std::complex<DP> LiebLin_ln_Overlap (Vect<std::complex<DP> > lstate_lambdaoc, DP lstate_lnnorm,
				       LiebLin_Bethe_State& rstate);

  // In src/LiebLin_Tgt0.cc:
  DP Entropy (LiebLin_Bethe_State& RefState);
  DP Canonical_Free_Energy (LiebLin_Bethe_State& RefState, DP kBT);
  LiebLin_Bethe_State Canonical_Saddle_Point_State (DP c_int, DP L, int N, DP kBT);
  LiebLin_Bethe_State Add_Particle_at_Center (const LiebLin_Bethe_State& RefState);
  LiebLin_Bethe_State Remove_Particle_at_Center (const LiebLin_Bethe_State& RefState);
  DP rho_of_lambdaoc_1 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta);
  DP rho_of_lambdaoc_2 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta);


} // namespace ABACUS

#endif