1#ifndef HILA_SUN_HEATBATH_H
2#define HILA_SUN_HEATBATH_H
27template <
typename T,
int N>
32 double xr1, xr2, xr3, xr4;
41 double pi2 = M_PI * 2.0;
49 action = U * staple.
dagger();
52 for (
int ina = 0; ina < N - 1; ina++)
53 for (
int inb = ina + 1; inb < N; inb++) {
61 v.d = action.
e(ina, ina).re + action.
e(inb, inb).re;
62 v.c = -(action.
e(ina, ina).im - action.
e(inb, inb).im);
63 v.a = -(action.
e(ina, inb).im + action.
e(inb, ina).im);
64 v.b = -(action.
e(ina, inb).re - action.
e(inb, ina).re);
105 d = -(xr2 + xr1 * xr3 * xr3) / al;
111 if ((1.00 - 0.5 * d) <= xr4 * xr4) {
114 for (k = 0; k < 40 && !test; k++) {
121 xr3 = cos(pi2 * xr3);
123 d = -(xr2 + xr1 * xr3 * xr3) / al;
124 if ((1.00 - 0.5 * d) > xr4 * xr4)
136 xl = exp((
double)(-2.0 * al));
140 for (k = 0; k < 40 && test == 0; k++) {
146 a0 = 1.00 + log((
double)r) / al;
147 if ((1.0 - a0 * a0) > xr2 * xr2)
169 r2 = 1.0 - a.d * a.d;
177 rho = r2 - a.c * a.c;
178 rho = sqrt(fabs(rho));
182 a.a = rho * cos((
double)xr2);
183 a.b = rho * sin((
double)xr2);
189 h2x2 = h.convert_to_2x2_matrix();
200 return (uhit / utry);
void mult_by_2x2_left(int p, int q, const Matrix_t< 2, 2, R, Mt > &M)
Multiply -matrix from the left by matrix defined by sub matrix.
Rtype dagger() const
Hermitian conjugate of matrix.
T e(const int i, const int j) const
Standard array indexing operation for matrices.
double random()
Real valued uniform random number generator.
double suN_heatbath(SU< N, T > &U, const SU< N, T > &staple, double beta)
heatbath
SU(N) Matrix definitions.