...

[phys/heatbath.git] / obs_phip2_hist.hpp

diff --git a/obs_phip2_hist.hpp b/obs_phip2_hist.hpp
index 437338338d582db34e6019e213783748ee5aa1d9..350c213ca1ec8835382bf36c5878b8d1d0a0f87b 100644 (file)
--- a/obs_phip2_hist.hpp
+++ b/obs_phip2_hist.hpp
@@ -8,6 +8,7 @@
 #include "latlib/obstat.hpp"
 
 #include <iostream>
+#include <stdlib.h>
 
 using namespace std;
 
@@ -15,57 +16,67 @@ class obs_phip2_hist : public o815::obs {
 
 public:
   struct obsmem {
-    double *phip2;
+    double phip2[4];
   };
-  obs_phip2_hist(o815 *_O815);
+  obs_phip2_hist(o815 *_O815, const int& avgmany=2);
   
 private:
   void _start() {}
   void _meas(bool loadedobs, const int& nthmeas);
   void _finish() {}
-  obsmem OM;
+  obsmem *OM;
   sim *Sim;
+  double avgphip2[4];
+  static string intToStr(const int& number);
+  int avgmany;
 };
 
-obs_phip2_hist::obs_phip2_hist(o815 *_O815) : o815::obs("phip2", 
-                                                       _O815->paraQ->getParaNames() + "nthstep:kpt:kpx:phip2", _O815, 
-                                                       sizeof(double)*_O815->comargs.lsize[0]*_O815->comargs.lsize[1], 
-                                                       "_hist") {
-  OM.phip2 = (double*)obsMem;
+string obs_phip2_hist::intToStr(const int& number)
+{
+  stringstream ss;
+  ss << number;
+  return ss.str();
+}
+
+obs_phip2_hist::obs_phip2_hist(o815 *_O815, const int& _avgmany) : o815::obs("phip2", 
+                                                                            _O815->paraQ->getParaNames() + "nthstep:imode:kt:kx:phip2absdiffequi", _O815, 
+                                                                            sizeof(obsmem), 
+                                                                            "_hist_avg" + intToStr(_avgmany) ) {
+  OM = (obsmem*)obsMem;
   Sim = (sim*)O815->Sim;
+  avgmany = _avgmany;
 }
 
 void obs_phip2_hist::_meas(bool loadedobs, const int& nthmeas) {
-  if (!loadedobs) {
-    for (int kpt = 0; kpt < O815->comargs.lsize[0]; kpt++)
-      for (int kpx = 0; kpx < O815->comargs.lsize[1]; kpx++) {
-       double ppt = 2.*M_PI/O815->comargs.lsize[0] * kpt;
-       double ppx = 2.*M_PI/O815->comargs.lsize[1] * kpx;
+  for (int kpimode=0; kpimode<4; kpimode++) {
+    const double ppt = 2.*M_PI/O815->comargs.lsize[0] * ( 0 + int( kpimode/4. * O815->comargs.lsize[0] ) ) ;
+    const double ppx = 2.*M_PI/O815->comargs.lsize[1] * ( 0 + int( kpimode/4. * O815->comargs.lsize[1] ) );
 
-       OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] = 0;
-       
-       for (int ixt = 0; ixt < O815->comargs.lsize[0]; ixt++)
-         for (int ixx = 0; ixx < O815->comargs.lsize[1]; ixx++)
-           for (int ixpt = 0; ixpt < O815->comargs.lsize[0]; ixpt++)
-             for (int ixpx = 0; ixpx < O815->comargs.lsize[1]; ixpx++)
-               OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] += real( conj( Sim->conf[ ixt*O815->comargs.lsize[1] + ixx ].phi )
-                                                                       * Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi
-                                                                       * exp ( _i_*(double)ppx*(double)(ixx-ixpx) + _i_*(double)ppt*(double)(ixt-ixpt) ) );
-       
-       OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] /= Sim->LSIZE2;
-      }
-  }
-  
-  for (int kpt = 0; kpt <= O815->comargs.lsize[0]; kpt++)
-    for (int kpx = 0; kpx <= O815->comargs.lsize[1]; kpx++) {
-      double ppt = 2.*M_PI/O815->comargs.lsize[0] * kpt;
-      double ppx = 2.*M_PI/O815->comargs.lsize[1] * kpx;
+    if (!loadedobs) {
+      complex<double> phitildepp = 0;
+      
+      for (int ixpt = 0; ixpt < O815->comargs.lsize[0]; ixpt++)
+       for (int ixpx = 0; ixpx < O815->comargs.lsize[1]; ixpx++)
+         phitildepp += Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi 
+           * exp ( - _i_*(double)ppx*(double)ixpx - _i_*(double)ppt*(double)ixpt );
+      
+      OM->phip2[ kpimode ] = norm( phitildepp) / Sim->LSIZE2;
+    }
+
+    if ( nthmeas%avgmany == 0 )
+      avgphip2[kpimode] = 0;
+    
+    avgphip2[kpimode] += OM->phip2[ kpimode ];
 
+    if ( (nthmeas+1)%avgmany == 0 ) {
       *out << O815->paraQ->getParaVals();
       *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
-      *out << "\t" << kpt << "\t" << kpx;
-      *out << "\t" << OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) << endl;
+      *out << "\t" << kpimode;
+      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[0] );
+      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[1] );
+      *out << "\t" << abs( avgphip2[kpimode] / (double)avgmany - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
     }
+  }
 };
 
 #endif