From 40af29b16e2d1054cb9e0433eee5a7f45fe811a5 Mon Sep 17 00:00:00 2001
From: Alexander Schmidt <alex@treefish.org>
Date: Thu, 28 Nov 2013 14:31:16 +0100
Subject: [PATCH 1/1] Averaging phip2 over 2 successive measurements.

---
 obs_phip2_hist.hpp | 41 ++++++++++++++++++++++++++++++-----------
 1 file changed, 30 insertions(+), 11 deletions(-)

diff --git a/obs_phip2_hist.hpp b/obs_phip2_hist.hpp
index 0f2ef32..350c213 100644
--- 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;
 
@@ -17,7 +18,7 @@ public:
   struct obsmem {
     double phip2[4];
   };
-  obs_phip2_hist(o815 *_O815);
+  obs_phip2_hist(o815 *_O815, const int& avgmany=2);
   
 private:
   void _start() {}
@@ -25,14 +26,25 @@ private:
   void _finish() {}
   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:imode:kt:kx:phip2absdiffequi", _O815, 
-							sizeof(obsmem), 
-							"_hist") {
+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) {
@@ -50,13 +62,20 @@ void obs_phip2_hist::_meas(bool loadedobs, const int& nthmeas) {
       
       OM->phip2[ kpimode ] = norm( phitildepp) / Sim->LSIZE2;
     }
+
+    if ( nthmeas%avgmany == 0 )
+      avgphip2[kpimode] = 0;
     
-    *out << O815->paraQ->getParaVals();
-    *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
-    *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( OM->phip2[ kpimode ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
+    avgphip2[kpimode] += OM->phip2[ kpimode ];
+
+    if ( (nthmeas+1)%avgmany == 0 ) {
+      *out << O815->paraQ->getParaVals();
+      *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
+      *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;
+    }
   }
 };
 
-- 
2.39.5