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

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

#define SQRT32 1.224744871391589
#define S0 0.7131271672829936
#define S1 9.63504163050617
#define S2 -11.750233478730895
#define S3 4.797012730243139


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 = velocity, v_T = thermal velocity
 */
{
	double e,v1,v2,x,y;
	do {
		do {
			do { 
				do v1=rand01(); while (v1==0.0); // should not be zero
				v2=2.0*rand01()-1.0;
			} while (v1*v1+v2*v2 > 1.0); 
			y=v2/v1; //y = tan(pi * U).
			x=S0*y+SQRT32; 
		} while (x <= 0.0); //Accept only x>0
		y = x*x;
		e = exp(-y)*y*x*(S1 + S2*x + S3*y); //Ratio of prob. fn. to comparison fn.		
	} while (rand01() > e); //Reject on basis of a second uniform deviate. 
	return x;
}

double Bounded_Maxwell_Velocity_Deviate(double x1, double x2)
/*
 * As above but x confined to x1<x<x2
 */
{
	double e,v1,v2,x,y;
	do {
		do {
			do { 
				do v1=rand01(); while (v1==0.0); // should not be zero
				v2=2.0*rand01()-1.0;
			} while (v1*v1+v2*v2 > 1.0); 
			y=v2/v1; //y = tan(pi * U).
			x=S0*y+SQRT32; 
		} while (x <= x1 || x >= x2); // Accept only within [x1, x2)
		y = x*x;
		e = exp(-y)*y*x*(S1 + S2*x + S3*y); //Ratio of prob. fn. to comparison fn.		
	} while (rand01() > e); //Reject on basis of a second uniform deviate. 
	return x;
}


#define NDEF 1000000
#define VMEAN 1.329340388179137
#define V2MEAN 2.0

int main(int argc, char* argv[])
{
	double s1 = 0, s2 = 0, v;
	int i,N;

	N = argc>1 ? atoi(argv[1]) : NDEF;

	srand( (unsigned int)time( NULL ) );

	printf("Testing the Maxwell random deviate generator.\n");
	printf("%d random deviates generated\n\n",N);
	
	for(i=1; i<=N; ++i)
	{
		v = Maxwell_Velocity_Deviate();
		s1 += v;
		s2 += v*v;
	}

	printf("Average velocity <v>.\n");
	printf("  Expected:   %g\n",VMEAN);
	printf("  Calculated: %g\n",s1/N);
	printf("  Rel. Deviation: %g\n\n",fabs(s1/N/VMEAN-1));

	printf("Average velocity squared <v^2>.\n");
	printf("  Expected:   %g\n",V2MEAN);
	printf("  Calculated: %g\n",s2/N);
	printf("  Rel. Deviation: %g\n",fabs(s2/N/V2MEAN-1));

	return 0;
}