#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#define rand01() ( ((double)rand())/((double)RAND_MAX+1) )

double Maxwell_Velocity_Deviate()
/*
 * Returns a random deviate x, 0<x<inf,
 * distributed as the maxwell flux distribution 2*exp(-x^2)*x^3, 
 * where x = v/v_T, v the velocity, v_T the thermal velocity
 */
{
	double u1,u2;
	do u1 = rand01(); while (u1==0.0);
	do u2 = rand01(); while (u2==0.0);
	return sqrt(-log(u1*u2));
}

/* default max number of iterations */
#define MAXN_DEF 1000000
/* limits for uniform velocity generation */
#define VMIN 0.1
#define VMAX 5.0

int main(int argc, char* argv[])
{
	double s1,s2,v,v2,w;
	int i, j;
	clock_t start;

	/* get user input */
	int MAXN = argc>1 ? atoi(argv[1]) : MAXN_DEF;

	printf(
		"Testing the Maxwell random deviate generator.\n\n"
		"The value of <v^2> = 2 will be calculated by 2 methods:\n"
		"  1. Generating random v distributed uniformly in %g < v <%g and\n"
		"     multiplying by the distribution weight w = 2*exp(-v^2)*v^3\n"
		"  2. Generating random v distributed according to 2*exp(-v^2)*v^3\n\n",
		VMIN,VMAX);

	/* seed the random generator */
	srand( (unsigned int)time( NULL ) );

	// 1. std method
	printf("1st method\n\n");
	printf("N          <v^2>      delta<v^2>/<v^2> time(ms)\n");
	s1 = s2 = 0.0;
	j = 100;
	start = clock();
	for(i=1; i<=MAXN; ++i)
	{
		/* generate a random velocity uniformly in [VMIN,VMAX) and calc v^2 */
		v = VMIN + rand01()*(VMAX-VMIN);  
		v2 = v*v;
		/* distribution weight */
		w = 2*exp(-v2)*v2*v;
		/* multiply by v^2 to calculate <v^2> */
		w *= v2;
		/* add w, w^2 to sums */
		s1 += w;
		w *= w;
		s2 += w;

		/* print iteration no, <v^2>, std. deviation and time */
		if (i==j)
		{
			j *= 10;
			printf("%10d %10.6f %16.3e %10.3f\n",i,
				(VMAX-VMIN)*s1/i,
				sqrt(s2/s1/s1-1./i),
				1.*(clock() - start) / CLOCKS_PER_SEC * 1000
				);
		}
	}

	printf("\n");


	// 2.  method
	printf("2nd method\n\n");
	printf("N          <v^2>      delta<v^2>/<v^2> time(ms)\n");
	s1 = s2 = 0.0;
	j = 100;
	start = clock();
	for(i=1; i<=MAXN; ++i)
	{
		/* generate a random velocity with the Maxwell flux distribution */
		v = Maxwell_Velocity_Deviate();  
		/* calc v^2 and add to the sums. weight = 1 */
		v2 = v*v;
		s1 += v2;
		s2 += v2*v2;

		/* print iteration no, <v^2>, std. deviation and time */
		if (i==j)
		{
			j *= 10;
			printf("%10d %10.6f %16.3e %10.3f\n",i,
				s1/i,
				sqrt(s2/s1/s1-1./i),
				1.*(clock() - start) / CLOCKS_PER_SEC * 1000
				);
		}
	}

	return 0;
}