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

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

Copyright (c) J.-S. Caux

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

File:  ABACUS_Usage_Example_Heis.cc

Purpose:  illustrates basic use of ABACUS for spin chains.

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

#include "ABACUS.h"

using namespace std;
using namespace JSC;


int main()
{
  clock_t StartTime = clock();

  // Refer to include/ABACUS_Heis.h for all class definitions
  // and to src/HEIS/ for the actual implementations.

  // Set basic system parameters:
  DP Delta = 1.0; // anisotropy
  int N = 64; // chain length
  int M = 32; // number of downturned spins; must be <= N/2 (on or above equator)

  // Define the chain:  J, Delta, h, Nsites
  Heis_Chain chain(1.0, Delta, 0.0, N);

  // The Heis_Chain class so constructed contains information about all types of strings:
  cout << "Delta = " << Delta << endl << "Possible string lengths and parities: ";
  for (int j = 0; j < chain.Nstrings; ++j)
    cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t";
  cout << endl;


  // Constructing a Bethe state: start with the ground state for simplicity.
  // Define the base:  chain, Mdown
  // A base by definition represents a partition of the M down spins into different strings.
  Heis_Base gbase(chain, M); // This puts all the M down spins into one-strings.

  // Define the ground state
  //XXZ_Bethe_State gstate(chain, gbase);  // Use this constructor for 0 < Delta < 1
  XXX_Bethe_State gstate(chain, gbase); // Use this constructor for Delta == 1
  //XXZ_gpd_Bethe_State gstate(chain, gbase); // Use this constructor for Delta > 1.
  // Anisotropies Delta < 0 are not separately implemented.

  // The assignment operator is overloaded for Bethe states:
  XXX_Bethe_State gstate2 = gstate; // copies all data of gstate into new object gstate2


  // Compute everything about the ground state
  gstate.Compute_All(true);

  // Output information about the state
  cout << gstate << endl;


  // Now define some excited state.

  // First method:
  // Start with defining a base using a base constructor with a vector of
  //  numbers of rapidities of each type as argument:
  // First define Nrapidities vector:
  Vect<int> Nrapidities(0, chain.Nstrings);
  // Assigning the numbers that follow, you have to ensure (yourself!) that
  //  the \sum_j M_j n_j = M (where n_j is the length of type j).
  Nrapidities[0] = M-2; // number of one-strings
  Nrapidities[1] = 1; // one two-string
  // Define the base:
  Heis_Base ebase(chain, Nrapidities);
  // Once the base is defined, the limiting quantum numbers are automatically computed:
  cout << "ebase defined, data is (Mdown, Nrap, Nraptot, Ix2_infty, Ix2_min, Ix2_max, baselabel):"
       << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl
       << ebase.Ix2_infty << endl
       << ebase.Ix2_min << endl << ebase.Ix2_max << endl << ebase.baselabel << endl;

  // An excited state can then be defined using this new base:
  XXX_Bethe_State estate(chain, ebase);
  // Individual quantum numbers can be manipulated: this will NOT update the
  //  state label or verify range of Ix2
  estate.Ix2[0][0] = M+1;
  estate.Compute_All(true);
  cout << endl << "estate: " << estate << endl;


  // Second method of constructing a base:
  Heis_Base ebase2(chain, "31");  // This is another constructor for all rapidities in one-strings
  //Heis_Base ebase(chain, "29x1y1"); // one two-string (XXX)


  // Construct a new Bethe state
  XXX_Bethe_State estateref (chain, ebase2); // defaults to the lowest-energy state for this base
  XXX_Bethe_State estate2 (chain, ebase2); // yet another state

  // Setting a state to a given label:
  // The base of estate must coincide with the base in label. Label is relative to estateref.Ix2.
  estate2.Set_to_Label ("31_1_nh", estateref.Ix2);
  estate2.Compute_All(true);
  cout << "estate2: " << estate2 << endl;

  // Energy and momentum of the states:
  cout << "Energy and momentum of states:" << endl;
  cout << "gstate.E = " << gstate.E << "\testate.E = " << estate.E << endl;
  cout << "Excitation energy = estate.E - gstate.E = " << estate.E - gstate.E << endl;
  cout << "Excitation momentum (mod 2 pi) = estate.K - gstate.K = " << estate.K - gstate.K << endl;


  // Computing matrix elements:
  // CAREFUL: magnetizations must be consistent with operator (error is flagged).
  cout << "Matrix elements: " << endl;
  // The logarithm of the matrix elements are computed as complex numbers:
  cout << "ln_Sz (gstate, estate) = " << ln_Sz_ME (gstate, estate) << endl;

  // The other possible matrix element call is for the Smin matrix element,
  cout << "ln_Smin (gstate, estate2) = " << ln_Smin_ME (gstate, estate2) << endl;


  clock_t StopTime = clock();

  //cout << "Total time: " << (StopTime - StartTime)/*/CLOCKS_PER_SEC*/ << " hundreths of a second."
  cout << "Total time: " << DP(StopTime - StartTime)/CLOCKS_PER_SEC << " seconds."
       << endl;

  return(0);
}