]> git.treefish.org Git - phys/heatbath.git/blob - sim-nor.hpp
2f5d22c50611f62c18463213220394b2f878608e
[phys/heatbath.git] / sim-nor.hpp
1 #ifndef SIM_HPP
2 #define SIM_HPP
3
4 #include <gsl/gsl_rng.h>
5 #include <complex>
6 #include <math.h>
7
8 #include "latlib/neigh.h"
9
10 #define EPSILONPHI 0.5
11
12 class sim : public o815::sim {
13 public:
14   struct siteconf {
15     complex<double> phi;
16   };
17   sim(o815 *_O815);
18   siteconf* conf;
19   unsigned int LSIZE2;
20
21 private:
22   void _makeSweep();
23   void _newParas();
24   gsl_rng* rangsl;
25   neigh *nb;
26   int updatePhi (const int& x);
27   double M;
28   double SofPhi(const int& x, const complex<double>& phix);
29   double rhoPhi(const int& x, const complex<double>& phixCandidate);
30 };
31
32
33
34 sim::sim(o815 *_O815) : o815::sim( _O815, 
35                                    sizeof(siteconf)*
36                                    (_O815->comargs.lsize[0]*_O815->comargs.lsize[1]) ) {
37   conf = (siteconf*)confMem;
38
39   rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
40   gsl_rng_set(rangsl, time(NULL));
41
42   LSIZE2 = _O815->comargs.lsize[0] * _O815->comargs.lsize[1];
43
44   nb = new neigh(2, _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
45 }
46
47 void sim::_makeSweep() {  
48   for (int ichecker=0; ichecker<2; ichecker++)
49     for (int it=0; it<O815->comargs.lsize[0]; it++)
50       for (int iy=(it+ichecker)%2; iy<O815->comargs.lsize[1]; iy+=2)
51         updatePhi( it*O815->comargs.lsize[1] + iy );
52 }
53
54 void sim::_newParas() {
55   M = pow( (*O815->paraQ)["mass"], 2 ) + 4;
56   *log << "SIM: Resetting fields." << endl << flush;
57   memset(conf, 0, sizeof(siteconf)*LSIZE2);
58 }
59
60 int sim::updatePhi(const int& x) {
61   complex<double> phixCandidate = conf[x].phi +
62     complex<double> ( 2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ), 
63                       2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ) );
64   
65   if ( gsl_rng_uniform(rangsl) < rhoPhi(x, phixCandidate) ) {
66     conf[x].phi = phixCandidate;
67     return 1;
68   }
69   
70   return 0;
71 }
72
73 double sim::rhoPhi(const int& x, const complex<double>& phixCandidate) {
74   return exp( SofPhi(x, conf[x].phi) - SofPhi(x, phixCandidate) );
75 }
76
77 double sim::SofPhi(const int& x, const complex<double>& phix) {
78   double sofphi = M * norm(phix);
79   
80   for (int nu=0; nu<4; nu++)
81     sofphi -= 2 * real( conj(phix) * conf[ (*nb)[x*4+nu] ].phi );
82   
83   return sofphi;
84 }
85
86 #endif