92 lines
3.5 KiB
C++
92 lines
3.5 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: LiebLin_ln_Overlap.cc
|
|
|
|
Purpose: Calculates Nikita Slavnov's determinant, case RHS is Bethe
|
|
and LHS is arbitrary
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
complex<DP> LiebLin_ln_Overlap (Vect<DP> lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate)
|
|
{
|
|
Vect<complex<DP> > lstate_lambdaoc_CX(rstate.N);
|
|
for (int i = 0; i < rstate.N; ++i) lstate_lambdaoc_CX[i] = complex<DP>(lstate_lambdaoc[i]);
|
|
return(LiebLin_ln_Overlap(lstate_lambdaoc_CX, lstate_lnnorm, rstate));
|
|
}
|
|
|
|
complex<DP> LiebLin_ln_Overlap (Vect<complex<DP> > lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate)
|
|
{
|
|
// Computes the log of the overlap between the left state and the Bethe rstate
|
|
|
|
// IMPORTANT ASSUMPTIONS:
|
|
// The length over which the overlap is calculated is that in rstate
|
|
|
|
if (lstate_lambdaoc.size() != rstate.N) {
|
|
cout << "Caution: calling overlap of states with difference charges. Are you sure that's what you want ?" << endl;
|
|
return(0.0);
|
|
}
|
|
|
|
// \mu are rstate rapidities, \lambda are from lstate (not BE)
|
|
|
|
Vect_CX ln_prod_plus(0.0, rstate.N); // contains \prod_a (\mu_a - lambdaoc_k + II)
|
|
Vect_CX ln_prod_minus(0.0, rstate.N); // contains \prod_a (\mu_a - \lambdaoc_k - II)
|
|
Vect_CX ln_prod_ll_plus (0.0, rstate.N); // contains \prod_{a \neq k} (\lambdaoc_a - \lambdaoc_k + II)
|
|
Vect_CX ln_prod_ll_minus (0.0, rstate.N); // contains \prod_{a \neq k} (\lambdaoc_a - \lambdaoc_k - II)
|
|
|
|
for (int j = 0; j < rstate.N; ++j)
|
|
for (int k = 0; k < rstate.N; ++k) {
|
|
ln_prod_plus[k] += log(rstate.lambdaoc[j] - lstate_lambdaoc[k] + II);
|
|
ln_prod_minus[k] += log(rstate.lambdaoc[j] - lstate_lambdaoc[k] - II);
|
|
if (j != k) {
|
|
ln_prod_ll_plus[k] += log(lstate_lambdaoc[j] - lstate_lambdaoc[k] + II);
|
|
ln_prod_ll_minus[k] += log(lstate_lambdaoc[j] - lstate_lambdaoc[k] - II);
|
|
}
|
|
}
|
|
|
|
// Build the matrix:
|
|
SQMat_CX Omega (0.0, rstate.N);
|
|
|
|
for (int j = 0; j < rstate.N; ++j)
|
|
for (int k = 0; k < rstate.N; ++k)
|
|
Omega[j][k] = exp(-II * rstate.cxL * lstate_lambdaoc[k] - ln_prod_plus[k] + ln_prod_minus[k])
|
|
* (II/((rstate.lambdaoc[j] - lstate_lambdaoc[k])*(rstate.lambdaoc[j] - lstate_lambdaoc[k] - II)))
|
|
- (II/((rstate.lambdaoc[j] - lstate_lambdaoc[k])*(rstate.lambdaoc[j] - lstate_lambdaoc[k] + II)));
|
|
|
|
complex<DP> lndetOmega = lndet_LU_CX_dstry(Omega);
|
|
|
|
// Prefactors:
|
|
complex<DP> ln_prod_d_mu = II * 0.5 * rstate.cxL * rstate.lambdaoc.sum();
|
|
complex<DP> ln_prod_d_lambdaoc = II * 0.5 * rstate.cxL * lstate_lambdaoc.sum();
|
|
|
|
complex<DP> ln_prod_mu = 0.0;
|
|
complex<DP> ln_prod_lambdaoc = 0.0;
|
|
complex<DP> ln_prod_plusminus = 0.0;
|
|
for (int j = 0; j < rstate.N - 1; ++j)
|
|
for (int k = j+1; k < rstate.N; ++k) {
|
|
ln_prod_mu += log(rstate.lambdaoc[k] - rstate.lambdaoc[j]);
|
|
ln_prod_lambdaoc += log(lstate_lambdaoc[k] - lstate_lambdaoc[j]);
|
|
}
|
|
|
|
for (int j = 0; j < rstate.N; ++j)
|
|
for (int k = 0; k < rstate.N; ++k)
|
|
ln_prod_plusminus += log((rstate.lambdaoc[j] - lstate_lambdaoc[k] + II));
|
|
|
|
return(ln_prod_d_mu + ln_prod_d_lambdaoc - ln_prod_mu - ln_prod_lambdaoc + ln_prod_plusminus
|
|
+ lndetOmega - 0.5 * (lstate_lnnorm + rstate.lnnorm));
|
|
}
|
|
|
|
} // namespace ABACUS
|