#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>
#include <vector>

using namespace std;

typedef double REAL;
typedef unsigned int uint;
#define TAB "\t"

int main()
{
	const REAL beta0 = 1.0050525387; //Beta with which calculations were done.
	const REAL rwradius = 0.1; //How 'long' to reweight from referencevalue.
	const int betacount = 101; //How many different betas to use.
	const uint skipcount = 1000; //How many measurements to skip.
	
	vector<REAL> measurements;	//Array where to store data. Data is read to memory for performance reasons.

	measurements.reserve(2*(200000-skipcount));
	//Reserve memory beforehand for the data for better performance.
	//In this case there is expected to be 200000 measurements in the file(|M| and E on every line).

	const REAL V = 64.0*64.0;
	ifstream fin("../atBetacL64.txt");
	//Open data file.
	//Note: The file is expected to contain nothing put the two data columns, |M| and E in this order.

	REAL temp;
	for(int i = 0; i<2*skipcount; i++) fin >> temp;
	//Skip measurements (number is given by skipcount).

	copy(istream_iterator<REAL>(fin), istream_iterator<REAL>(), back_inserter(measurements));
	//Read data from file.

	fin.close();

	cout << "Number of measurements: " << measurements.size()/2 << endl;

	for(int i = 0; i<betacount; i++)
	{
		const REAL beta = beta0 - rwradius + i * 2.0 * rwradius / (betacount - 1.0);
		//Calculate current beta value. Values are taken from [beta0 - rwradius, beta0 + rwradius] with
		//even step size.

		const REAL betaDiff = beta - beta0;
		REAL  denominator = 0, Mexp = 0, M2exp = 0;

		//Go through the data one row at a time.
		for(vector<REAL>::const_iterator iter = measurements.begin(); iter != measurements.end(); iter++)
		{
			const REAL M = *iter; iter++;
			const REAL E = *iter;
			const REAL eval = exp(-betaDiff * E);
			Mexp += M * eval;
			M2exp += M * M * eval;
			denominator += eval;
		}
		Mexp /= denominator;
		M2exp /= denominator;
		cout << beta << TAB << Mexp << TAB << M2exp << TAB
			 << (1.0/V) * (M2exp - Mexp*Mexp) << endl;
	}

	return 0;
}