for (int nu=0; nu<4; nu++)
V += conf[ (*nb)[x*4+nu] ].phi;
-
+
const double V2diff = pow(real(V), 2) - pow(imag(V), 2);
const double Vprod = real(V)*imag(V);
conf[x].phi = complex<double> ( + real(conf[x].phi) * V2diff + 2 * imag(conf[x].phi) * Vprod,
- imag(conf[x].phi) * V2diff + 2 * real(conf[x].phi) * Vprod ) / norm(V);
-
+
const double r = exp( - M * norm(conf[x].phi)
- 1./M * norm(V)
+ 2 * real( conf[x].phi * conj(V) ) );
conf[x].phi = sqrt(std::log( 1./(1-r) )) / sqrt(M)
- * polar( 1.0, arg(conf[x].phi) )
+ * polar( 1.0, gsl_rng_uniform(rangsl) * 2*M_PI )
+ V / M;
}