Only computing IR and UV mode.

[phys/heatbath.git] / sim-nor.hpp

diff --git a/sim-nor.hpp b/sim-nor.hpp
index aa305ff71abcaaec1e0e1a7f976e671ffe13d74d..b67f8a783b5d2c4e1b7f33107ddfcadcbb245d8e 100644 (file)
--- a/sim-nor.hpp
+++ b/sim-nor.hpp
@@ -7,6 +7,8 @@
 
 #include "latlib/neigh.h"
 
+#define EPSILONPHI 0.5
+
 class sim : public o815::sim {
 public:
   struct siteconf {
@@ -15,16 +17,17 @@ public:
   sim(o815 *_O815);
   siteconf* conf;
   unsigned int LSIZE2;
+  double M;
+  double m;
+  neigh *nb;
 
 private:
   void _makeSweep();
   void _newParas();
-
   gsl_rng* rangsl;
-
-  neigh *nb;
-
-  void updatePhi (const int& x);
+  int updatePhi (const int& x);
+  double SofPhi(const int& x, const complex<double>& phix);
+  double rhoPhi(const int& x, const complex<double>& phixCandidate);
 };
 
 
@@ -42,17 +45,46 @@ sim::sim(o815 *_O815) : o815::sim( _O815,
   nb = new neigh(2, _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
 }
 
-void sim::updatePhi (const int& x) 
-{
-  const double r = gsl_rng_uniform(rangsl);
-  const double theta = gsl_rng_uniform(rangsl) * 2*M_PI;
-}
-
 void sim::_makeSweep() {  
+  for (int ichecker=0; ichecker<2; ichecker++)
+    for (int it=0; it<O815->comargs.lsize[0]; it++)
+      for (int iy=(it+ichecker)%2; iy<O815->comargs.lsize[1]; iy+=2)
+       updatePhi( it*O815->comargs.lsize[1] + iy );
 }
 
 void sim::_newParas() {
-  /* reset variables */
+  m = (*O815->paraQ)["mass"];
+  M = pow( (*O815->paraQ)["mass"], 2 ) + 4;
+  *log << "SIM: Resetting fields." << endl << flush;
+
+  for (int ix=0; ix<LSIZE2; ix++)
+    conf[ix].phi = 1;
+}
+
+int sim::updatePhi(const int& x) {
+  complex<double> phixCandidate = conf[x].phi +
+    complex<double> ( 2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ), 
+                     2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ) );
+  
+  if ( gsl_rng_uniform(rangsl) < rhoPhi(x, phixCandidate) ) {
+    conf[x].phi = phixCandidate;
+    return 1;
+  }
+  
+  return 0;
+}
+
+double sim::rhoPhi(const int& x, const complex<double>& phixCandidate) {
+  return exp( SofPhi(x, conf[x].phi) - SofPhi(x, phixCandidate) );
+}
+
+double sim::SofPhi(const int& x, const complex<double>& phix) {
+  double sofphi = M * norm(phix);
+  
+  for (int nu=0; nu<4; nu++)
+    sofphi -= 2 * real( conj(phix) * conf[ (*nb)[x*4+nu] ].phi );
+  
+  return sofphi;
 }
 
 #endif