/* This is a program for implementing a site-bond percolation model that I developed. Instead of
assigning a fixed bonding arrangement with a fixed coordination number, this model allows for
random bond assignment with a fixed coordination number within a cubic lattice. It also allows
for the attachment of multiple molecules at a single bond.  The program executes without any 
command line inputs, and queries for relevant information.  It outputs some data to the screen 
and all data to datafile. 
Author: Riyad Chetram Raghu Date:November 15, 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 <fstream.h>
#include <time.h>
#include <Classes.cpp>
#include <percfunc.cpp>

void main()
{
	cout << "please input the lattice parameter L:\n";  // this prompts the user for inputs
	cin >> L;
	cout << "please input the site occupancy probability:\n";
	float psite;
	cin >> psite;
	cout << "please input the bonding probability:\n";
	cin >> pbond;
	cout << "please input the number of PNIPAAm segments per molecule:\n"; 
	short CN;
	cin >> CN;
	cout << "please wait...\n" << flush; 
	sample.InitLattice();  // this initializes the lattice
	srand((unsigned)time(NULL));
	long rcluscount=0,ncluscount=0;
	float montecarlo;

	Site *thissite;
	for (k=0;k<L;k++)
	{
		for (j=0;j<L;j++)
		{
			for (i=0;i<L;i++) // these three nested forloops visit every site in the lattice
			{
				thissite=sample.getsite(i,j,k);
				(*thissite).initsite();  // initialize the site at i,j,k
				montecarlo=RandomReal(0,1); // produce a random real number between 0 and 1;
				if (montecarlo<=psite)  // In this case a molecule occupies the site
				{
					rcluscount++;
					if(realstart==&nill) // In this case no molecules have been assigned yet
					{  
						current=new Cluster;
						realstart=current;
						(*current).setnext(pnill);
						(*current).setprev(pnill);
					}
					else
					{
						current=AddCluster(current);
					};
					(*thissite).cluster=current;
					(*thissite).bondcount=CN;
					InitCluster(rcluscount,i,j,k,current);
				}
				else // In this case the site is empty
				{
					ncluscount++;
					if(emptystart==&nill) // in this case no sites have been assigned yet
					{
						ncurrent=new Cluster;
						emptystart=ncurrent;
						(*ncurrent).setnext(pnill);
						(*ncurrent).setprev(pnill);
					}
					else
					{
						ncurrent=AddCluster(ncurrent);
					};
					(*thissite).cluster=ncurrent;
					(*thissite).bondcount=-1;
					InitCluster(ncluscount,i,j,k,ncurrent);
				};
			};
		};
	};
	cout << "The lattice has now been filled with molecules.\n" << flush;


	// this is the segment of the algorithm that forms bonds
	current=realstart; // assign the current pointer to look at the first cluster
	while (current!=&nill)
	{
		i=(*current).x[0]; // i,j and k are set to look at the root site of the cluster
		j=(*current).y[0];
		k=(*current).z[0];
		BondAlg((*current).label); // the root site of the current cluster is examined to form new bonds
		current=(*current).getnext();  // the algorithm moves on to the next cluster
	};
	cout << "The molecule containing clusters have now been formed.\n" << flush;  // debuging statement
// the code works to this point

	// this is the segment of the algorithim that forms empty space clusters
	ncurrent=emptystart; // assign the current pointer to look at the first empty site
	while (ncurrent!=&nill)
	{
		i=(*ncurrent).x[0]; // i,j and k are set to look at the root site of the cluster
		j=(*ncurrent).y[0];
		k=(*ncurrent).z[0];
		NullAlg((*ncurrent).label); // All of the empty sites connected to the root site are added to the current cluster
		ncurrent=(*ncurrent).getnext(); // the algorithm moves on to the next cluster
	};

	// the lattice has now been filled and following code deals with statistics collection
	cout << "The empty site clusters have now been formed.\n";
	ofstream fout("data.txt",ios::noreplace);
	if (!fout)
	{  cout << "This would rewrite the existing data.txt file\n";
	   return;
	};

	int *distribution=new int[(L*L*L)]; // this array will store the distribution of cluster sizes
	for (long a=0;a<L*L*L;a++) distribution[a]=0;
	double bigcluster=0;
	double gelfraction=0; // this variable will store the sum of the members in the percolating clusters and eventually the gel fraction
	double radgyr=0; // this variable will store the radius of gyration;
	double memsum=0,memsqsum=0,corrlength=0,mem=0;
	
	// This section of code gets the statistics of the molecule containing clusters
	current=realstart;
	while (current!=&nill)
	{
		// this if statement checks for percolation
		if(((*current).z[1]==0&&(*current).z[2]==L-1)||((*current).y[1]==0&&(*current).y[2]==L-1)||((*current).x[1]==0&&(*current).x[2]==L-1)) 
		{
			gelfraction=gelfraction +(*current).members;  // this adds up all of the molecules in a percolating cluster
			fout << 1 << "\t";
			cout << "Percolation has occured.\n";
		}
		else fout << 0 << "\t";
		if (bigcluster<(*current).members) bigcluster=(*current).members;
		mem=(*current).members;
		memsum=mem+memsum;
		memsqsum=mem*mem+memsqsum;
		fout << (*current).members << "\t"; // this outputs the number of molecules in the cluster
		// this is the radius of gyration for the current cluster
		radgyr=(mem*(*current).rsquared -(*current).xsum*(*current).xsum-
			(*current).ysum*(*current).ysum-(*current).zsum*(*current).zsum)/(mem*mem);
		corrlength=corrlength+mem*mem*radgyr;
		fout << radgyr << "\n"; // this outputs the radius of gyration of the cluster
		distribution[(*current).members-1]=distribution[(*current).members-1]+1;
		current=(*current).getnext(); // this looks at the next cluster
	};
	gelfraction=gelfraction/memsum; // the gel fraction of the simulation
	corrlength=2*corrlength/memsqsum; // the correlation length for the sumulation
	cout << "The gel fraction is:" << gelfraction << "\n";
	cout << "The erodable fraction is:" << 1-gelfraction << "\n";
	cout << "The correlation length is:" << corrlength << "\n";
	cout << "The weight average cluster size is:" << memsqsum/memsum << "\n";
	cout << "The weight average cluster size without the largest cluster is:" << (memsqsum-bigcluster*bigcluster)/(memsum-bigcluster) << "\n";
	fout << "\n" << gelfraction << "\n";
	fout << 1-gelfraction << "\n";
	fout << corrlength << "\n";
	fout << memsqsum/memsum << "\n";
	fout << (memsqsum-bigcluster*bigcluster)/(memsum-bigcluster) << "\n\n";
	for (a=0;a<L*L*L;a++)
	{
		fout << a+1 << "\t" << distribution[a] << "\n"; // this writes the size distribution of the clusters to file
		distribution[a]=0; // this resets the distribution for the next section
	};
	fout <<"\n";

	// This section of code gets the statistics of the empty clusters
	memsum=0;
	memsqsum=0;
	corrlength=0;
	gelfraction=0;
	ncurrent=emptystart;
	while (ncurrent!=&nill)
	{
		// this if statement checks for percolation
		if(((*ncurrent).z[1]==0&&(*ncurrent).z[2]==L-1)||((*ncurrent).y[1]==0&&(*ncurrent).y[2]==L-1)||((*ncurrent).x[1]==0&&(*ncurrent).x[2]==L-1)) 
		{
			fout << 1 << "\t";
			cout << "Empty site percolation has occured.\n";
		}
		else fout << 0 << "\t";
		mem=(*ncurrent).members;
		memsum=mem+memsum;
		memsqsum=mem*mem+memsqsum;
		fout << (*ncurrent).members << "\t"; // this outputs the number of molecules in the cluster
		// this is the radius of gyration for the current cluster
		radgyr=(mem*(*ncurrent).rsquared -(*ncurrent).xsum*(*ncurrent).xsum-
			(*ncurrent).ysum*(*ncurrent).ysum-(*ncurrent).zsum*(*ncurrent).zsum)
			/(mem*mem);
		corrlength=corrlength+mem*mem*radgyr;
		fout << radgyr << "\n"; // this outputs the radius of gyration of the cluster
		distribution[(*ncurrent).members-1]=distribution[(*ncurrent).members-1]+1;
		ncurrent=(*ncurrent).getnext(); // this looks at the next cluster
	};
	corrlength=2*corrlength/memsqsum; // the correlation length for the sumulation
	cout << "The correlation length of the empty clusters is:" << corrlength << "\n";
	cout << "The weight average empty cluster size is:" << memsqsum/memsum << "\n";
	fout << "\n" << corrlength << "\n";
	fout << memsqsum/memsum << "\n\n";
	for (a=0;a<L*L*L;a++)
	{
		fout << a+1 << "\t" << distribution[a] << "\n"; // this writes the size distribution of the clusters to file
	};
	fout.close();
	cout << "Hit enter to exit.\n" << flush;
	cin.ignore(1,'\n');
	cin.get();
}