//Tool for using the Flyvbjerg-Petersen blocking method to calculate the error on the sampling averages. This file contains code for calculating the error on the mean potential energy. By reading a plot of this code's output, one should observe a plateau, which is the estimate of the error on the mean potential energy
#include <vector.h>
#include <math.h>
#include <stdio.h>

int main() { 
  
   FILE * energyData = fopen( "Energy.dat", "r" );

   vector<double> kineticEnergies;
   vector<double> potentialEnergies;
   
   double kineticEnergy = 0.0;
   double potentialEnergy = 0.0;
   double totalEnergy = 0.0; 

   int count = 0;

   while ( !feof( energyData ) ) {
      fscanf( energyData, "%lf %lf %lf", &kineticEnergy, &potentialEnergy, &totalEnergy );  
      if ( !feof( energyData ) ) {
	kineticEnergies.push_back( kineticEnergy );
	potentialEnergies.push_back( potentialEnergy );
	count++;
      }
   }

   fclose( energyData );

   printf ( "%d\n", count );

   count /= 2;

   FILE * constantFile = fopen( "ErrorEstimate.dat", "w" );

   for ( int i=1; i<count; i*= 2 ) {

     const int itmp = i;
     double block[count/itmp];
     for ( int j=0; j<count/i; j++ ) {
       block[j] = 0.0;
     }

     double average = 0.0;
     double meanSquared = 0.0;

     for ( int j=0; j<count/i; j++ ) {
       
       for ( int k=count+j*i; k<count+j*i+i; k++ ) {
	 block[j] += potentialEnergies[k];
       }
       block[j] /= ((double)(i));
       average += block[j];
       meanSquared += block[j]*block[j];
     }

     average /= ((double)(count/i));
     meanSquared /= ((double)(count/i));
     
     double stdDev = 0.0;

     for ( int j=0; j<count/i; j++ ) {
       stdDev += (block[j]-average)*(block[j]-average);
     }
     stdDev /= ((double)(count/i));

     printf ( "i = %d, nblocks = %d, mean squared = %lf, standard dev = %lf, average = %lf, <x^2>-<x>^2 = %lf, the constant = %20.14e\n", i, count/i, meanSquared, stdDev, average, meanSquared-average*average, (meanSquared-average*average)/((double)((count/i)-1)) );
     fprintf ( constantFile, "%e  %e\n", sqrt((meanSquared-average*average)/((double)((count/i)-1))),  sqrt((meanSquared-average*average)/((double)((count/i)-1)))*(1./sqrt(2.*((double)((count/i)-1)))) );
     
   }

   fclose ( constantFile );


}