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

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

Copyright (c) J.-S. Caux.

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

File:  ABACUS_XXX_VOA.h

Purpose:  Declares classes for XXX in zero field: Vertex Operator Approach

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

#ifndef ABACUS_XXX_VOA_H
#define ABACUS_XXX_VOA_H

#include "ABACUS.h"

const DP default_req_prec = 1.0e-14;
const int default_max_rec = 10;

namespace ABACUS {

  inline DP Integrand_11 (Vect_DP args)
  {
    // args[0] corresponds to t, args[1] to rho
    return((exp(args[0]) * sinh(args[0]) * cos(4.0 * args[0] * args[1])
	    /(2.0 * pow(cosh(args[0]), 2.0)) + 2.0 * pow(sin(2.0 * args[0] * args[1]), 2.0))/args[0]);
  }

  inline DP Integrand_12 (Vect_DP args)
  {
    DP expm2t = exp(-2.0*args[0]);
    return(cos(4.0 * args[0] * args[1]) * expm2t * (3.0 + expm2t)/ (args[0] * (1.0 + expm2t) * (1.0 + expm2t)));
  }

  inline DP Integrand_2 (Vect_DP args)
  {
    DP answer = 0.0;
    if (args[0] < 1.0) answer = exp(args[0]) * pow(sin(2.0 * args[0] * args[1]), 2.0)
			 /(args[0] * sinh(args[0]) * pow(cosh(args[0]), 2.0));
    else if (args[0] >= 1.0) {
      DP expm2t = exp(-2.0 * args[0]);
      answer = 8.0 * expm2t * pow(sin(2.0 * args[0] * args[1]), 2.0)
	/(args[0] * (1.0 - expm2t) * (1.0 + expm2t) * (1.0 + expm2t));
    }
    return(answer);
  }

  inline DP Integrand_A (Vect_DP args)
  {
    // This kernel is for -ln | A_-(i\pi/2) | function
    return(exp(args[0]) * pow(sinh(args[0]/2.0), 2.0)/(args[0] * sinh(2.0*args[0]) * cosh(args[0])));
  }

  DP I_integral (DP rho, DP req_prec);


  /*********************  TWO SPINONS ********************/

  DP SF_2p (DP k, DP omega, I_table Itable);
  DP SF_2p (Vect_DP args, I_table Itable);
  DP SF_2p_alt (Vect_DP args, I_table Itable);
  DP SF_2p_w (Vect_DP args, I_table Itable);
  DP SF_2p_w_alt (Vect_DP args, I_table Itable);
  DP SF_2p_intw (Vect_DP args, I_table Itable);
  DP SF_2p_intw_alt (Vect_DP args, I_table Itable);
  DP SF_2p_check_sumrule (DP req_prec, int max_rec, I_table Itable);
  DP SF_2p_check_sumrule_alt (DP req_prec, int max_rec, I_table Itable);
  DP Fixed_k_sumrule_w (DP k);
  DP Fixed_k_sumrule_omega (DP k);
  DP SF_2p_check_fixed_k_sumrule (DP k, DP req_prec, int max_rec, I_table Itable);
  DP SF_2p_check_fixed_k_sumrule_alt (DP k, DP req_prec, int max_rec, I_table Itable);
  DP SF_2p_check_fixed_k_sumrule_opt (DP k, DP req_prec, int Npts, I_table Itable);


  /********************** FOUR SPINONS **********************/

  DP Sum_norm_gl (Vect_DP rho, DP req_prec);
  DP Compute_C4 (DP req_prec);
  DP SF_contrib (Vect_DP p, DP req_prec, I_table Itable);
  DP J_fn (Vect_DP p, DP req_prec, I_table Itable);
  inline DP Jacobian_p3p4_KW (DP k, DP w, DP K, DP W)
  {
    return(1.0/sqrt(pow(twoPI * sin(0.5 * (k - K)), 2.0) - (w-W)*(w-W)));
  }
  bool Set_p_given_kwKW (DP k, DP w, DP K, DP W, Vect_DP& p);
  inline DP wmin_4p (DP k)
  {
    return(PI * fabs(sin(k)));
  }
  inline DP wmax_4p (DP k)
  {
    return(2.0 * PI * sqrt(2.0 * (1.0 + fabs(cos(0.5*k)))));
  }
  inline  DP Wmin (DP k, DP w, DP K)
  {
    return(ABACUS::max(1.0e-15, ABACUS::max(fabs(PI * sin(K)), w - twoPI * sin(0.5 * (fabs(k-K))))));
  }
  inline  DP Wmax (DP k, DP w, DP K)
  {
    return(ABACUS::min(twoPI * sin(0.5 * K), w - fabs(PI * sin(k - K))));
  }
  DP G_fn (Vect_DP args_to_G, I_table Itable);
  DP G1_fn (Vect_DP args_to_G, I_table Itable);
  DP G2_fn (Vect_DP args_to_G, I_table Itable);
  DP G1_fn_mid (Vect_DP args_to_G, I_table Itable);
  DP G2_fn_mid (Vect_DP args_to_G, I_table Itable);
  DP G_fn_alt (Vect_DP args_to_G, I_table Itable);
  DP H_fn (Vect_DP args_to_H, I_table Itable);
  DP H2_fn (Vect_DP args_to_H, I_table Itable);
  DP H_fn_mid (Vect_DP args_to_H, I_table Itable);
  DP H_fn_alt (Vect_DP args_to_H, I_table Itable);
  DP SF_4p_kwKW (Vect_DP args, I_table Itable);
  DP SF_4p_kwKW_alpha (Vect_DP args, I_table Itable);
  DP SF_4p_kwKW_alpha_opt (Vect_DP args, I_table Itable);

  // Interface to used version:
  DP SF_4p_rec (DP k, DP omega, DP req_prec, int max_rec, I_table Itable);
  DP SF_4p (DP k, DP omega, I_table Itable);
  DP SF_4p_opt (DP k, DP omega, DP req_prec, int Npts_K, int Npts_W, I_table Itable);
  void Translate_raw_4p_data (DP k, int max_rec_w, const char* SFraw_Cstr, const char* SF_Cstr, const char* SFsrc_Cstr, I_table Itable);
  DP SF_4p_rec (DP k, DP req_prec, int max_rec_w, int max_rec, I_table Itable);
  Integral_result SF_4p_opt (DP k, DP req_prec, int Npts, I_table Itable);
  Integral_result SF_4p_opt (DP k, DP req_prec, int Npts_w, int Npts_KW, I_table Itable);
  Integral_result SF_4p_opt (DP k, DP req_prec, int Npts_w, int Npts_K, int Npts_W, I_table Itable);

  //******************************** Functions to produce files similar to ABACUS **********************************
  void Produce_SF_2p_file (int N, int Nomega, DP omegamax, I_table Itable);
  void Produce_SF_4p_file (int N, int Nomega, DP omegamax, DP req_prec, int max_rec, I_table Itable);
  void Produce_SF_4p_file (int N, int Nomega, DP omegamax, DP req_prec, int Npts_K, int Npts_W, I_table Itable);

  //******************************** New version, full k and omega integral in one go ******************************
  DP Direct_J_integral (int Npts_p, DP req_prec, I_table Itable);
  DP Direct_J_integral_bin (int Npts_p, int Npts_o, DP req_prec, I_table Itable);
  void Smoothen_raw_SF_4p (int Npts_p, int Npts_o, DP req_prec, DP width);

} // namespace ABACUS

#endif