From: Alexander Schmidt Date: Tue, 26 Nov 2013 10:17:29 +0000 (+0100) Subject: Fixed wrong 1mr+ update. X-Git-Url: http://git.treefish.org/~alex/phys/heatbath.git/commitdiff_plain/2845de363d157f4e7b3ad37c4d4b13a8d77bf449?hp=274a66cdb22fd36ea8cbbc373783d848f1ba253e Fixed wrong 1mr+ update. --- diff --git a/sim-1mr+.hpp b/sim-1mr+.hpp index 74fc18e..e46a9ea 100644 --- a/sim-1mr+.hpp +++ b/sim-1mr+.hpp @@ -24,9 +24,6 @@ private: void _newParas(); gsl_rng* rangsl; void updatePhi (const int& x); - void updatePhi_angle_relax (const int& x, const complex& V); - void updatePhi_angle_rand (const int& x, const complex& V); - void updatePhi_r (const int& x, const complex& V); }; @@ -44,30 +41,6 @@ sim::sim(o815 *_O815) : o815::sim( _O815, nb = new neigh(2, _O815->comargs.lsize[0], _O815->comargs.lsize[1]); } -void sim::updatePhi_r (const int& x, const complex& 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) ) - + V / M; -} - -void sim::updatePhi_angle_relax (const int& x, const complex& V) -{ - const double V2diff = pow(real(V), 2) - pow(imag(V), 2); - const double Vprod = real(V)*imag(V); - conf[x].phi = complex ( + real(conf[x].phi) * V2diff + 2 * imag(conf[x].phi) * Vprod, - - imag(conf[x].phi) * V2diff + 2 * real(conf[x].phi) * Vprod ) / norm(V); -} - -void sim::updatePhi_angle_rand (const int& x, const complex& V) -{ - conf[x].phi = polar( abs(conf[x].phi), gsl_rng_uniform(rangsl) * 2*M_PI ); -} - void sim::updatePhi (const int& x) { complex V=0; @@ -75,12 +48,18 @@ void sim::updatePhi (const int& x) for (int nu=0; nu<4; nu++) V += conf[ (*nb)[x*4+nu] ].phi; - updatePhi_r(x,V); + const double V2diff = pow(real(V), 2) - pow(imag(V), 2); + const double Vprod = real(V)*imag(V); + conf[x].phi = complex ( + real(conf[x].phi) * V2diff + 2 * imag(conf[x].phi) * Vprod, + - imag(conf[x].phi) * V2diff + 2 * real(conf[x].phi) * Vprod ) / norm(V); - if ( gsl_rng_uniform(rangsl) < 0.5 ) - updatePhi_angle_relax(x,V); - else - updatePhi_angle_rand(x,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, gsl_rng_uniform(rangsl) * 2*M_PI ) + + V / M; } void sim::_makeSweep() {