// Program file name: GMC.java                                        
// Monte Carlo simulation of grain growth
// SI units are used throughout the code

import java.lang.*;
import java.io.*;
import java.util.Random;
public class GMC {
  
  static Random r = new Random();
  static double xc=0.0070001,yc=0.0050001;   //location of the torch
  static int nx=60,ny=60;
  static double x[]=new double[nx];
  static double y[]=new double[ny];
  static double t[][]=new double[nx][ny];
  static double tmcs[][]=new double[nx][ny];
  static double hx,hy,xmax=0.01,ymax=0.01;    //grid size and dimension of plate
  static double diameter;
  static double vel=0.005;
  
  public static void main(String argv[]) {
  
	hx=xmax/nx;
	hy=ymax/ny;
    for(int i=0;i<nx;i++)
	    for(int j=0;j<ny;j++)
		{
		    x[i]=i*hx;
		    y[j]=j*hy;
		 }
  
    callroutine();
  }
   
//initializes orientation, calls temperature, calls tmc and calls potts 
public static void callroutine(){

		int o[][]=new int[nx][ny];
		double avg=0,maxt=0.;
			
		for(int i=0;i<nx;i++)
			for(int j=0;j<ny;j++)
			{
				o[i][j]=Math.abs(r.nextInt()%32);
			}
		temperature(t);
		displayD(t,"temperature.dat");
		tmcs=tmc(t);										 //monte-carlo time at each location
	
		for(int i=0;i<nx;i++)
			for(int j=0;j<ny;j++){
	    			avg+=tmcs[i][j];
	    			maxt=(tmcs[i][j]>maxt)?tmcs[i][j]:maxt;
	    			}
	   avg/=nx*ny;
		
		System.out.println("\nmax and average tmcs value: " + maxt + "  " +avg);
		displayD(tmcs,"tmcs.dat");
	
		for(int k=0;k<=20;k++)
			potts(o,t,tmcs);
				
		diameter=grainsize(o); 
		display(o);
   }

public static void temperature(double t[][])    //quasi-steady thermal profile
    {
		double rad, al;
		double q0=500,d=0.002;			//net power, thickness of plate
		double rho=7800,cp=600,k=45;    //density,specific heat cpacity, thermal conductivity
		double t0=298.0;				//initial temperature
		al=k/(rho*cp);					//thermal diffusivity
		
		for(int i=0;i<nx;i++)
			for(int j=0;j<ny;j++)
			{
				rad=Math.sqrt((x[i]-xc)*(x[i]-xc)+(y[j]-yc)*(y[j]-yc));
				t[i][j]=t0+((q0/d)/(2*3.14159*k))*Math.exp(-0.5*vel*(x[i]-xc)/al)*bessk0(0.5*vel*rad/al);
			}
    }

public static double[][] tmc(double t[][])      // monte-carlo time at each location (i,j)
	{
		double tmcs[][]=new double[nx][ny];
		double expt[][]=new double[nx][ny];
		double time[]=new double[nx];
		double K1=1.01;
		double A=1;
		double Vm=7.11e-6;
		double h=6.62e-34;		    //Planck's constant
		double lam=hx;
		double L0=hx;
		double Sf=9.48;
		double R=8.314;
		double Q=9.34e4;
		double n1=0.42;
		double gam=1.77;
		double N=1.53e19;
		double Na=6.022e23;			//Avogadro's number
		double dt,C1;
		C1=(4*gam*A*N*Vm*Vm)/(Na*Na*h)*Math.exp(Sf/R);
		
		for(int i=0;i<nx;i++)
			time[i]=x[i]/vel;
		
		for(int j=0;j<ny;j++){
				for(int i=nx-1;i>=0;i--){
				dt=(i<nx-1)?time[i+1]-time[i]:time[i];
				expt[i][j]=(i<nx-1)?expt[i+1][j]+Math.exp(-Q/(R*t[i][j]))*dt:Math.exp(-Q/(R*t[i][j]))*dt;
			
				}
		}
		
		for(int j=0;j<ny;j++)
				for(int i=0;i<nx;i++)
					tmcs[i][j]=Math.pow((Math.sqrt(C1*expt[i][j]+L0*L0)/(K1*lam)-1/K1),(1/n1));
		return tmcs;
    }
    
public static void potts(int o[][],double t[][],double tmcs[][])  //update orientation using potts model for energy reduction
{
		int ii,jj,ic,E1,E2;
		double prob,maxt;
		int oc[][]=new int[nx][ny];
		
		maxt=tmcs[0][0];
		
		for(int i=0;i<nx;i++)
			for(int j=0;j<ny;j++){
	    			oc[i][j]=o[i][j];
	    			maxt=(tmcs[i][j]>maxt)?tmcs[i][j]:maxt;
	    			}
		
		for(int k=0;k<50000000;k++){
			
	  		ii=Math.abs(r.nextInt()%nx);
			jj=Math.abs(r.nextInt()%ny);
			prob=tmcs[ii][jj]/maxt;
		
			if(prob>Math.abs(r.nextDouble())){
				E1=energy(ii,jj,o);
				ic=Math.abs(r.nextInt()%32);
				oc[ii][jj]=ic;
				E2=energy(ii,jj,oc);	
				if (E2<E1) o[ii][jj]=ic;
				}	
		}
	}
	
public static int energy(int ii,int jj,int o[][])    //energy of location (ii,jj)
{
		int i,j,E=0;

		i=(ii+1)%nx;  	j=jj;
		E-=(o[i][j]==o[ii][jj])?0:-1;
		
		i=(ii+1)%nx;  	j=(jj+1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1;
		
		i=(ii+1)%nx;  	j=(jj+ny-1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1; 
		
		i=ii;  	j=(jj+1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1;

		i=ii;  	j=(jj+ny-1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1;

		i=(ii+nx-1)%nx;  	j=jj;
		E-=(o[i][j]==o[ii][jj])?0:-1;
		
		i=(ii+nx-1)%nx;  	j=(jj+1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1;

		i=(ii+nx-1)%nx;  	j=(jj+ny-1)%ny;
		E-=(o[i][j]==o[ii][jj])?0:-1;
		
		return E;
	
	}
	
public static double grainsize(int o[][])     //measurement of grain size distribution and average grain size 
  {
		int dist[]=new int[10];
		double avg,size;
		int sum=0;
		int num=0;
		
		for(int i=0;i<10;i++)
			dist[i]=0;
			
		for(int i=0;i<nx;i+=5)					//measurement of average grain size using mean lineal intercept method
			for(int j=0;j<ny;j++){
				if (o[i][j]!=o[i][(j+ny-1)%ny]) // periodic boundary
					num+=1;
					}
		avg=nx*ymax/num;						//average grain size
		
		for(int i=0;i<nx;i++)					// grain size distribution
			for(int j=0;j<ny;j++){
				if (o[i][j]==o[i][(j+ny-1)%ny])   
					sum+=1;
				else
					{
					size=sum*hy/avg;
					sum=0;	
					for(int k=0;k<10;k++)
					{
						if (size<((k<9)?(k+1):1000)*0.5 && size>=k*0.5)
							dist[k]+=1;
					}
					}
				}
		System.out.println("\nGrain size distribution :\n  ");	
		for(int i=0;i<10;i++)
			System.out.println(i*0.5 + " - " + ((i<9)?(i+1):1000)*0.5 + "   " +dist[i] );	
		System.out.println("\nAverage grain size :  " + avg );				
		return avg ;	
  }

public static void displayD(double tmcs[][],String filename)   //write temperature/tMCS to file for viewing using software Tecplot 
{
	FileOutputStream out; // declare a file output object
    PrintStream p; // declare a print stream object
			 try
				{
                       out = new FileOutputStream(filename);
                        p = new PrintStream( out );
						p.println ("Time distribution");
						p.println("VARIABLES = 'X', 'Y', 't'");
						p.println("ZONE I="+ nx+ "J=" +ny+ "F=POINT");
						for(int j=0;j<ny;j++)
							for(int i=0;i<nx;i++)
	    							p.println(x[i]+"  " +y[j]+"   "+tmcs[i][j]);
						p.close();
                }
                catch (Exception e)
                {
                        System.err.println ("Error writing to file");
                }
}
  
 
public static void display(int o[][])    //write orientation to file for viewing using software Tecplot
  {
   			    
	   FileOutputStream out; 
       PrintStream p; 
       try
                {
                
                        out = new FileOutputStream("grain.dat");
                        p = new PrintStream( out );
						p.println ("Grain orientation");
						p.println("VARIABLES = 'X', 'Y', 'O'");
						p.println("ZONE I="+ nx+ "J=" +ny+ "F=POINT");
						for(int j=0;j<ny;j++)
							for(int i=0;i<nx;i++)
	    							p.println(x[i]+"  " +y[j]+"   "+o[i][j]);
						p.close();
				}
				catch (Exception e)
				{
						System.err.println ("Error writing to file");
				}
		
  }
  
public static double bessi0(double x)    // Returns the modified Bessel function I0(x) for any real x.
	{
			double ax,ans;
			double y; 
			if ((ax=Math.abs(x)) < 3.75) 
				{
					y=x/3.75;
					y*=y;
					ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
					+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
				} 
			else{
					y=3.75/ax;
					ans=(Math.exp(ax)/Math.sqrt(ax))*(0.39894228+y*(0.1328592e-1
					+y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
					+y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
					+y*0.392377e-2))))))));
				}
			return ans;
	}

public static double bessk0(double x)   // Returns the modified Bessel function K0(x) for positive real x.
	{
		  double y,ans; 
		  if (x <= 2.0)
			 {
				y=x*x/4.0;
				ans=(-Math.log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
				+y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
				+y*(0.10750e-3+y*0.74e-5))))));
			 } 
		  else
			 {
				y=2.0/x;
				ans=(Math.exp(-x)/Math.sqrt(x))*(1.25331414+y*(-0.7832358e-1
				+y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
				+y*(-0.251540e-2+y*0.53208e-3))))));
			  }
		  return ans;
	} 
	
}