/*
 ============================================================================
 Name        : main.c
 Author      : Racz Eva
 Version     : 1.0
 Copyright   : BME-TTK FI SzamSzim gyak segedanyag
 Description : Fail of naive Metropolis on Maxwell-Boltzmann in C, Ansi-style
 ============================================================================
 */

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(void) {
	const double beta = 1.0;	// temperature
	const int N = 1E3;			// number of particles
	double energy[N], dE;		// energy
	int i,j,k;					// indices
	double dim = 2.0;			// dimension of the box
	int D = (int) floor(dim);	// dimension of the box

	gsl_rng* r;

	r = gsl_rng_alloc(gsl_rng_mt19937);
	gsl_rng_set(r, 645864125);

	/* Initial energy of the particles */
	double E = 0.0;
	for (i=0; i < N; i++) {
		energy[i] = 2.0*(0.5/beta)*gsl_ran_chisq(r, dim);
		E += energy[i];
	}

	double energy_scale = 0.5/beta/2.0;
	double E_uj;

	/* Do a million steps changing the _momenta_ */
	for (i=0; i<1E6; i++) {
		/* Chose a random particle */
		j = gsl_rng_uniform_int(r,N);
		/* and its new energy */
		E_uj = 2.0*(0.5/beta)*gsl_ran_chisq(r, dim);

		/*
		E_uj = 0.0;
		for(k = 0; k < D; ++k){
			E_uj += pow(sqrt(energy[j]/dim)+(2.0*gsl_rng_uniform(r)-1.0)*sqrt(energy_scale), 2.0);
		}
		dE = E_uj - energy[j];
		*/


		if ( gsl_rng_uniform(r) < exp(-beta*dE) ) {
			energy[j] += dE;
			E += dE;
	//      printf(" 1 %lf\n", dE);
		} else {
	//      printf(" 0 %lf\n", dE);
		}
		printf("%d\t%g\n", i, E);
	}

	/* Output */
	for(i=0; i < N; i++)
		fprintf(stderr, "%g\n", energy[i]);

	gsl_rng_free(r);
	return 0;
}