/* This file contains the necessary files for the application of the percolation algorithm, this includes look up tables, 
vertice and neighbour search functions, and the random number generator.
Author: Riyad Chetram Raghu		Date: November 16, 2001 */
// Copyright (c) 2002-2003 Riyad Chetram Raghu

// Version 1.0 of January 13, 2003.
 

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>

int i=0,j=0,k=0;
short vertices[8]={0,0,0,0,0,0,0,0};
short bonds[26]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
Cluster *currsite=&nill;
float pbond=0.9;
Lattice sample;
Cluster *current=pnill,*ncurrent=pnill,*realstart=pnill,*emptystart=pnill;

short bondltable[26][3]= {	{-1,+1,-1},  // a look up table to convert to find a neighbour's location form its bond code
							{0,+1,-1},
							{+1,+1,-1},
							{-1,0,-1},
							{0,0,-1},
							{+1,0,-1},
							{-1,-1,-1},
							{0,-1,-1},
							{+1,-1,-1},
							{-1,+1,0},
							{0,+1,0},
							{+1,+1,0},
							{-1,0,0},
							{+1,0,0},
							{-1,-1,0},
							{0,-1,0},
							{+1,-1,0},
							{-1,+1,+1},
							{0,+1,+1},
							{+1,+1,+1},
							{-1,0,+1},
							{0,0,+1},
							{+1,0,+1},
							{-1,-1,+1},
							{0,-1,+1},
							{+1,-1,+1}};

short vertltable[8][14]= {	{1,2,4,5,10,11,13,8,7,6,5,4,3,2}, // a look up table for finding the neighbours of a vertice
							{2,3,5,6,11,12,14,8,7,6,5,4,3,1}, // and which is the equivalent vertice for that neighbour
							{4,5,7,8,13,15,16,8,7,6,5,4,2,1},
							{5,6,8,9,14,16,17,8,7,6,5,3,2,1},
							{10,11,13,18,19,21,22,8,7,6,4,3,2,1},
							{11,12,14,19,20,22,23,8,7,5,4,3,2,1},
							{13,15,16,21,22,24,25,8,6,5,4,3,2,1},
							{14,16,17,22,23,25,26,7,6,5,4,3,2,1}};

float RandomReal(float minimum,float maximum)
{
	float a=rand(),temp=RAND_MAX;
	a=a/temp;
	a=a*(maximum-minimum)+minimum;
	return a;
}

int RandomInt(float minimum, float maximum)
{
	float a=rand(),temp=RAND_MAX;
	a=a/temp*(maximum-minimum)+minimum;
	int b=a;
	temp=b;
	if(a-temp>=0.5) b++;
	return b;
}

void getneighbour(short neighbour) // a function to change the cartesian coordinate globals to point at the neighbour site
{ 
	neighbour=neighbour-1;
	int a=0,b=0,c=0;
	a=bondltable[neighbour][0]+i;
	b=bondltable[neighbour][1]+j;
	c=bondltable[neighbour][2]+k;
	if (a>=0&&a<L) i=a;
	else // in this case looping is occuring
	{ 
		if (a<0) i=L-1; // this is looping from the x=0 to x=L-1
		else i=0;  //this is looping from the x=L-1 to x=0
	};
	if (b>=0&&b<L) j=b;
	else // in this case looping is occuring
	{
		if (b<0) j=L-1; // this is looping from the y=0 to y=L-1
		else j=0; //this is looping from the y=L-1 to y=0
	};
	if (c>=0&&c<L) k=c;
	else // in this case looping is occuring
	{
		if (c<0) k=L-1; // this is looping from the z=0 to z=L-1
		else k=0; //this is looping from the z=L-1 to z=0
	};
	return;
};

short getparent(short neighbour) //a function to find the parent from the perspective of the neighbour 
{
	return 27-neighbour;
};

short getvertex(short vertex,short neighbour)  //a function to find the vertex from the perspective of the neighbour
{
	vertex=vertex-1;
	for(short a=0;a<6;a++)
	{
		if(vertltable[vertex][a]==neighbour) return vertltable[vertex][a+7];
	};
	cout << "Error: This neighbour does not border on the vertex \n";
	return 0;
};

void Construct() //a function to construct the vertices table and the bond table
{	
	Site *rootsite;
	int l=i,m=j,n=k;
	for(short a=1;a<9;a++)
	{
		rootsite=sample.getsite(l,m,n);
		if ((*rootsite).checkvert(a)==1)  //in this case a bond has already been formed at this vertice
		{
			a=a-1;
			for(short b=0;b<7;b++)
			{
				if ((*rootsite).checkbond(vertltable[a][b])) //if the bond has already been formed only add
				{												 //bonds that have been marked as bonded to the 
					vertices[a]=vertices[a]+1;							 //vertices array
					bonds[vertltable[a][b]-1]=1;	// set these molecules to a code of 1 for the bonds array
				};
			};
			a=a+1;
		}
		else //in this case no bonds have been formed
		{
			a=a-1;
			for (short b=0;b<7;b++)
			{
				getneighbour(vertltable[a][b]); // if the bond has not already been formed only add neighbours with
				rootsite=sample.getsite(i,j,k);
				if ((*rootsite).bondcount>0)// free bonds to vertices
				{	
					vertices[a]=a+1;
					if(abs(i-l)>1||abs(j-m)>1||abs(k-n)>1) bonds[vertltable[a][b]-1]=1; //looping has occured
					else if(bonds[vertltable[a][b]-1]!=1) bonds[vertltable[a][b]-1]=2; // if the bond code is not 1 change it to 2
				};
				i=l; //reset i,j and k because getneighbour changes them
				j=m;
				k=n;
			};
			a=a+1;
		};
	};
	for(a=1;a<8;a++) vertices[a]=vertices[a]+vertices[a-1]; // this loop sets up vertices for selection by random numbers
	return;
}

void formbond(short vertex, short num) // this is a function to from a bond at a vertice
{	Site *rootsite,*neighbour;
	rootsite=sample.getsite(i,j,k);
	(*rootsite).bondcount--; // this line decrement the number of bonding segments available
	if ((*rootsite).checkvert(vertex)) return; // if this bond has been bonded before the algorithm exits
	(*rootsite).changevert(vertex,1); // this line marks this vertex as being bonded
	int l=i,m=j,n=k; // this stores a backup copy of this sites x,y,z coordinates
	short b=0,c=0;
	float d=0;
	int vert[7]={0,0,0,0,0,0,0};
	vertex=vertex-1;
	for (short a=0;a<7;a++) // this loop sets up the vert array for random neighbour selection
	{
		getneighbour(vertltable[vertex][a]);  // this sets the i,j,k coordinates to point at the neighbour in vertltable
		neighbour=sample.getsite(i,j,k);
		if ((*neighbour).bondcount>0) b++;
		vert[a]=b;
		(*neighbour).changevert(vertltable[vertex][a+7],1); // this statement marks that an attempt has been made to bond this vertex
		i=l; // resets the i,j,k variables to the coordinates of the current site
		j=m;
		k=n;
	};

	if (b==0) // in this case there is nothing to bond to so the bondcount and the vertices arrays must be restored
	{
		(*rootsite).bondcount++; // this resets bondcount
		(*rootsite).changevert(vertex+1,0); // this statement and the next for loop the vertex boolean arrays are reset
		for (a=0;a<7;a++)
		{
			getneighbour(vertltable[vertex][a]);
			neighbour=sample.getsite(i,j,k);
			(*neighbour).changevert(vertltable[vertex][a+7],0);
			i=l;
			j=m;
			k=n;
		}
		short e; // this last code segment changes the vertices array so that this vertex is not accessed again		
		if (vertex==0) e=vertices[vertex];
		else e=vertices[vertex]-vertices[vertex-1];
		for (a=vertex;a<8;a++) vertices[a]=vertices[a]-e;
		return;
	};

	for (num;num>0;num--)
	{
		c=RandomInt(1,b);// produce a random integer c between 1 and b
		for(a=0;c>vert[a];a++);
		d=RandomReal(0,1);//produce a random real number between 0 and 1;
		if(d<=pbond) // a bond has been formed
		{
			c=vertltable[vertex][a];
			getneighbour(c);
			neighbour=sample.getsite(i,j,k);
			(*rootsite).changebond(c,1); // changes the current site's bond bitfield to indicate this neighbour has been bonded
			(*neighbour).bondcount--; // decreases the neighbours bond count
			c=getparent(c); 
			(*neighbour).changebond(c,1); //changes the neighbour site's bond bitfield to indicate this nieghbour has been bonded
			i=l; //resets i,j,k to the coordinates of the current site
			j=m;
			k=n; 
			b--; //decreases the number of bonds that have to be checked
			if (b==0) return;
			for (a=a;a<7;a++) vert[a]=vert[a]-1; //this loop prevents the same neighbour from being accessed again;
		};
	};
	return;
};

void BondFormer(long clusnum,bool isroot) // this algorithm actually forms all of the bonds for a site
{
	for (short a=0;a<26;a++) bonds[a]=0;  // reinitialize the bonds and vertices arrays
	for (short b=0;b<8;b++) vertices[b]=0;
	Site *rootsite=sample.getsite(i,j,k);
	currsite=(*rootsite).cluster;  //this gets the cluster of the current site 
	if((*currsite).label==clusnum&&isroot==false) return; // if the algorithm has already been run on this site exit
	if((*currsite).label!=clusnum) 
	{
		MergeCluster(currsite,current); //otherwise merge this site's cluster to the current cluster
		(*rootsite).cluster=current; //change the site's cluster pointer to point at the root of the site
	};
	if((*rootsite).bondcount<1) return; //if all the bonds have been assigned exit
	Construct();  //build the necessary arrays for bond formation

	while ((*rootsite).bondcount>0) // This loop assigns bonds to vertices
	{
		if (vertices[7]==0) (*rootsite).bondcount=0; // In this case there are no neighbours left to bond
		else
		{
			a=RandomInt(1,vertices[7]); // produce a random integer between 0 and the value of vertices 7
			for (b=0;vertices[b]<a;b++); // this finds the first value of vertices which is greater than or equal to a
			a=a-vertices[b-1]; // this is the diference between a and the last value of vertices which is less than a 
			formbond(b+1,a);  // This command forms the bonds
		};
	};
}

void BondAlg(long clusnum) //this is the bonding algorithm
{	
	Site *rootsite=sample.getsite(i,j,k);
	BondObj *currcall=new BondObj; //this is the first member of the BondObj linked list
	BondObj *lastcall=currcall; // lastcall holds the address of the last member of the list
	BondObj *nextcall=currcall;
	(*currcall).setnext(pneg); // the next and prev pointers of the fist member are set to &neg
	(*currcall).setprev(pneg);
	(*currcall).initBondObj(i,j,k); //the site is initiated;
	BondFormer(clusnum,true);  // the bonds for the site are formed

	for(short a=0;a<26;a++) if((*rootsite).checkbond(a+1)==true&&bonds[a]==2) (*currcall).setneighbour(a,true);
	// this for loop sets up the neighbour bitfield for the next step

	short b;
	while(nextcall!=&neg) //this while loop runs BondFormer on the other members of the cluster
	{
		i=(*currcall).getx(); // i,j and k are reset to the values of the current site
		j=(*currcall).gety();
		k=(*currcall).getz();
		//this for loop runs bondformer on all of the neighbours of the current site
		for (b=0;b<26;b++) if ((*currcall).getneighbour(b)==true)  // if the neighbour is bonded, BondFormer is run on it
		{
			getneighbour(b+1);  // the x,y and z coordinates of the neighbour are retrieved
			rootsite=sample.getsite(i,j,k);
			lastcall=AddBondObj(lastcall); // a new BondObj is added to the end of the linked list
			(*lastcall).initBondObj(i,j,k); // the BondObj is initialized
			BondFormer(clusnum,false);
			for(a=0;a<26;a++) if((*rootsite).checkbond(a+1)==true&&bonds[a]==2) (*lastcall).setneighbour(a,true); 
			// this for loop sets up the neighbour bifield for this site  
			i=(*currcall).getx(); // i,j and k are reset to the values of the current site
			j=(*currcall).gety();
			k=(*currcall).getz();
		};
		
		// the current BondObj is no longer needed so it can be deleted
		nextcall=(*currcall).getnext();
		DelBondObj(currcall);
		currcall=nextcall;		
	};
	return;
}

void NullFormer(long clusnum, bool isroot) // this is the algorithm that actually connects the empty sites
{	Site *neighbour, *rootsite=sample.getsite(i,j,k);
	for (short a=0;a<26;a++) bonds[a]=0;
	currsite=(*rootsite).cluster; //this gets the cluster of the current site
	if((*currsite).label==clusnum&&isroot==false) return; // if the algorithm has already been run on this site exit
	if((*currsite).label!=clusnum) 
	{
		MergeCluster(currsite,ncurrent); // otherwise merge this site's cluster to the current cluster
		(*rootsite).cluster=ncurrent;
	};
	int l=i,m=j,n=k; //This stores the coordinates of this site in active memory
	for (a=0;a<26;a++)
	{
		getneighbour(a+1); // this visits each neighbour
		neighbour=sample.getsite(i,j,k);
		if((*neighbour).bondcount==-1&&(*neighbour).checkbond(26-a)==false) // the site is empty and it hasn't been added yet
		{
			(*neighbour).changebond(26-a,true); // change the bond bitfields to record the bond
			(*rootsite).changebond(a+1,true);
			if (abs(i-l)<2&&abs(j-m)<2&&abs(k-n)<2)	bonds[a]=1; //note that this site has to added to the cluster
		};
		i=l;  // this resets i,j and k to the coordinates of this site 
		j=m;
		k=n;
	};
}

void NullAlg(long clusnum) //this is the bonding algorithm for empty sites
{
	BondObj *currcall=new BondObj; //this is the first member of the BondObj linked list
	BondObj *lastcall=currcall; // lastcall holds the address of the last member of the list
	BondObj *nextcall=currcall;
	(*currcall).setnext(pneg); // the next and prev pointers of the fist member are set to &neg
	(*currcall).setprev(pneg);
	(*currcall).initBondObj(i,j,k); //the site is initiated;
	NullFormer(clusnum,true);  // the empty sites surrounding the current site are found

	for(short a=0;a<26;a++) if(bonds[a]==1) (*currcall).setneighbour(a,true);
	// this for loop sets up the neighbour bitfield for the next step

	short b;
	while(nextcall!=&neg) //this while loop runs NullFormer on the other members of the cluster
	{
		i=(*currcall).getx(); // i,j and k are reset to the values of the current site
		j=(*currcall).gety();
		k=(*currcall).getz();
		//this for loop runs NullFormer on all of the neighbours of the current site
		for (b=0;b<26;b++) if ((*currcall).getneighbour(b)==true)  // if the neighbour is bonded, BondFormer is run on it
		{
			getneighbour(b+1);  // the x,y and z coordinates of the neighbour are retrieved
			lastcall=AddBondObj(lastcall); // a new BondObj is added to the end of the linked list
			(*lastcall).initBondObj(i,j,k); // the BondObj is initialized
			NullFormer(clusnum,false);
			for(a=0;a<26;a++) if(bonds[a]==1) (*lastcall).setneighbour(a,true); 
			// this for loop sets up the neighbour bifield for this site  
			i=(*currcall).getx(); // i,j and k are reset to the values of the current site
			j=(*currcall).gety();
			k=(*currcall).getz();
		};
		
		// the current BondObj is no longer needed so it can be deleted
		nextcall=(*currcall).getnext();
		DelBondObj(currcall);
		currcall=nextcall;		
	};
	return;
}