#include <stdio.h>
#include <string.h>
#include <math.h>
// Initiallization routine for random number
void   RandomInitialise(int,int);
// Call to generate a random number between 0  and 10
double RandomUniform(void);
// Call to generate a Gaussian random number
double RandomGaussian(double,double);
// Call to generate a random integer between the two arguments
int    RandomInt(int,int);
// call to generate a random double between the two arguments
double RandomDouble(double,double);
FILE *ofp;
int main()
{
	int i=0; int j=0; double k=0; double m=4; int N = 300; double c=0; int steps=10;  int a =0; int b =0; int x =0; int y =0;
	double Water_x[1000] = {0.0}; double Water_y[1000]={0.0}; double D = 0.0; double temp1_x =0.0; double temp1_y =0.0;
	double Radius= 1; double Inner_Radius = 0.01; double sigma = 3.165; double epsilon = 0.6502;
	double temp2_x = 0.0; double temp2_y = 0.0; double d_x =0.0; double d_y =0.0;
	double E =0; double Random_Energy = 0.0; double b0 = 1.0; int g=0; int h=0;

	ofp = fopen("tweak.txt","w");
	RandomInitialise(1802,9373);
	////Initialize a Lattice for all the molecules to begin from....

	for(i=0; i<N; i++)
	{
		Water_x[i] = c;
		Water_y[i] = m--;
		if(m==0)
		{
			m=4;
			c++;
		}
	}


	////////NEED RADIUS OF GYRATION!!!!!!!!!!!!!!!!!!!!
	for(k=0; k<N*N*N;k++)
	{
		//particle picked at random
		a = RandomInt(0,N);
		b = RandomInt(0,N);
		//Original Value stored in temporary variable.
		temp1_x = Water_x[a];
		temp1_y = Water_y[b];
		//position changed at random
		Water_x[a] = b0*RandomDouble(-2,2);
		Water_y[b] = b0*RandomDouble(-2,2);
		//	fprintf(ofp, "New position %f, %f\n", Water_x[a], Water_y[b]);

		//Distance calculation for neighboring particles.

		for(i=0; i<N; i++)
		{
			d_x = pow(((Water_x[a] + Radius) - Water_x[i]),2);
			d_y = pow(( (Water_y[b] + Radius) - Water_y[i]),2);
			D = sqrt(d_x + d_y);

			if(D < Inner_Radius)
			{
				Water_x[a] = temp1_x;
				Water_y[b] = temp1_y;
				break;
			}
			if(D<2*Radius)
			{
				E = 4*epsilon*(pow((sigma/D),12) - pow((sigma/D),6));
				Random_Energy = RandomInt(0,500000000000);
				if(Random_Energy<E)
				{
					Water_x[a] = temp1_x;
					Water_y[b] = temp1_y;
				}
			}
		}
	}
	k=0;
	for(i=0; i<N; i++)
	{
		fprintf(ofp,"%f\n", Water_x[i]);
		//fprintf(ofp,"\n");
		//fprintf(ofp,"%f\n", Water_y[i]);
	}
	fprintf(ofp,"\n");
	for(i=0; i<N; i++)
	{
		fprintf(ofp,"%f\n", Water_y[i]);
	}



}










#define FALSE 0
#define TRUE 1

/*

This Random Number Generator is based on the algorithm in a FORTRAN
version published by George Marsaglia and Arif Zaman, Florida State

University; ref.: see original comments below.

At the fhw (Fachhochschule Wiesbaden, W.Germany), Dept. of Computer
Science, we have written sources in further languages (C, Modula-2
Turbo-Pascal(3.0, 5.0), Basic and Ada) to get exactly the same test
results compared with the original FORTRAN version.
April 1989
Karl-L. Noell <[email protected]>
and  Helmut  Weber <[email protected]>

This random number generator originally appeared in "Toward a Universal
Random Number Generator" by George Marsaglia and Arif Zaman.
Florida State University Report: FSU-SCRI-87-50 (1987)
It was later modified by F. James and published in "A Review of Pseudo-
random Number Generators"

THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.

(However, a newly discovered technique can yield
a period of 64^600. But that is still in the development stage.)
It passes ALL of the tests for random number generators and has a period
of 2^144, is completely portable (gives bit identical results on all
machines with at least 24-bit mantissas in the floating point
representation).

The algorithm is a combination of a Fibonacci sequence (with lags of 97
and 33, and operation "subtraction plus one, modulo one") and an
"arithmetic sequence" (using subtraction).


Use IJ = 1802 & KL = 9373 to test the random number generator. The
subroutine RANMAR should be used to generate 20000 random numbers.
Then display the next six random numbers generated multiplied by 4096*4096
If the random number generator is working properly, the random numbers

should be:

6533892.0  14220222.0  7275067.0

6172232.0  8354498.0   64633180.0

*/


/* Globals */

double u[97],c,cd,cm;

int i97,j97;

int test = FALSE;


/*

This is the initialization routine for the random number generator.
NOTE: The seed variables can have values between:    0 <= IJ <= 31328
0 <= KL <= 30081

The random number sequences created by these two seeds are of sufficient
length to complete an entire calculation with. For example, if sveral
different groups are working on different parts of the same calculation,
each group could be assigned its own IJ seed. This would leave each group
with 30000 choices for the second seed. That is to say, this random
number generator can create 900 million different subsequences -- with
each subsequence having a length of approximately 64^30.

*/

void RandomInitialise(int ij,int kl)

{

	double s,t;

	int ii,i,j,k,l,jj,m;


	/*

	Handle the seed range errors

	First random number seed must be between 0 and 31328

	Second seed must have a value between 0 and 30081

	*/

	if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {

		ij = 1802;

		kl = 9373;

	}


	i = (ij / 177) % 177 + 2;

	j = (ij % 177)       + 2;

	k = (kl / 169) % 178 + 1;

	l = (kl % 169);


	for (ii=0; ii<97; ii++) {

		s = 0.0;

		t = 0.5;

		for (jj=0; jj<24; jj++) {

			m = (((i * j) % 179) * k) % 179;

			i = j;

			j = k;

			k = m;

			l = (53 * l + 1) % 169;

			if (((l * m % 64)) >= 32)

				s += t;

			t *= 0.5;

		}

		u[ii] = s;

	}


	c    = 362436.0 / 16777216.0;

	cd   = 7654321.0 / 16777216.0;

	cm   = 16777213.0 / 16777216.0;

	i97  = 97;

	j97  = 33;

	test = TRUE;

}


/*

This is the random number generator proposed by George Marsaglia in

Florida State University Report: FSU-SCRI-87-50

*/

double RandomUniform(void)

{

	double uni;


	/* Make sure the initialisation routine has been called */

	if (!test)

		RandomInitialise(1802,9373);


	uni = u[i97-1] - u[j97-1];

	if (uni <= 0.0)

		uni++;

	u[i97-1] = uni;

	i97--;

	if (i97 == 0)

		i97 = 97;

	j97--;

	if (j97 == 0)

		j97 = 97;

	c -= cd;

	if (c < 0.0)

		c += cm;

	uni -= c;

	if (uni < 0.0)

		uni++;


	return(uni);

}


/*

ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.

THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,

VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

The function returns a normally distributed pseudo-random number

with a given mean and standard devaiation.  Calls are made to a

function subprogram which must return independent random

numbers uniform in the interval (0,1).

The algorithm uses the ratio of uniforms method of A.J. Kinderman

and J.F. Monahan augmented with quadratic bounding curves.

*/

double RandomGaussian(double mean,double stddev)

{

	double  q,u,v,x,y;


	/*

	Generate P = (u,v) uniform in rect. enclosing acceptance region

	Make sure that any random numbers <= 0 are rejected, since

	gaussian() requires uniforms > 0, but RandomUniform() delivers >= 0.

	*/

	do {

		u = RandomUniform();

		v = RandomUniform();

		if (u <= 0.0 || v <= 0.0) {

			u = 1.0;

			v = 1.0;

		}

		v = 1.7156 * (v - 0.5);


		/*  Evaluate the quadratic form */

		x = u - 0.449871;

		y = fabs(v) + 0.386595;

		q = x * x + y * (0.19600 * y - 0.25472 * x);


		/* Accept P if inside inner ellipse */

		if (q < 0.27597)

			break;


		/*  Reject P if outside outer ellipse, or outside acceptance region */

	} while ((q > 0.27846) || (v * v > -4.0 * log(u) * u * u));


	/*  Return ratio of P's coordinates as the normal deviate */

	return (mean + stddev * v / u);

}


/*

Return random integer within a range, lower -> upper INCLUSIVE

*/

int RandomInt(int lower,int upper)
{
	return((int)(RandomUniform() * (upper - lower + 1)) + lower);
}


/*

Return random float within a range, lower -> upper

*/
double RandomDouble(double lower,double upper){return((upper - lower) * RandomUniform() + lower);}