// Monte Carlo simulation with the
// Metropolis scheme for spinodal decomposition

import java.lang.*;
import java.util.*;
import java.io.*;
import java.awt.*;  // uses the abstract windowing toolkit

public class Flip {
  
  // change the following globals to suit your fancy:
  static int maxSize = 400;      // width of window in pixels
  static int dim = 100;          // lattice size; can be any divisor of maxSize
  static int squareWidth = maxSize / dim;
  static Graphics g;  // the frame's graphics object, for drawing into
  static int frameTop, frameLeft;  // platform-dependent width of frame's top and left edges
  static Color upColor = new Color(00,160,60);
  static Color downColor = new Color(255,245,220);
 	
  static int seed;	
  static final int nsize = 400;
  static final int nskip = 100;
  static final int ntotal = nsize*nskip;
  static final int neq = 0;		//# of iterations for equilibriation
  static final int lx = dim;	//dimensions of lattice
  static final int ly = dim;	//dimensions of lattice
  static int iaccept = 0;	// counter for calculating the acceptance rate of attempted flips
  static int at;
  static double r1, r2;
  static double T1 = 1.5;	//temperature at which simulation is to be done
  static int a[][] = new int[lx][ly];	// array that stores the spins in the lattice
  static Random r = new Random();
  static double Teq = 1.0;	//for equilibriation at temperature Teq before final simulations
  static double T;
  static double J = 1.0; 
  static double fn = 0.5; //fraction of negative spin sites
  static double fp = 1 - fn; //fraction of positive spin sites
  
  public static void main(String argv[]) {
  	// initialize window and graphics:
        Frame gWin = new Frame("Flip Model");
        gWin.setLocation(50,50);  // screen coordinates of top left corner
        gWin.setResizable(false);
        gWin.setVisible(true);  // show it!
        Insets theInsets = gWin.getInsets();
        gWin.setSize(maxSize+theInsets.left+theInsets.right,maxSize+theInsets.top+theInsets.bottom);
        frameTop = theInsets.top;
        frameLeft = theInsets.left;
        long resumeTime = System.currentTimeMillis() + 1000;
        do {} while (System.currentTimeMillis() < resumeTime);  // wait for frame to get initialized
        g = gWin.getGraphics();  // stick graphics object into a global so we can draw from anywhere!

  	
  	FileOutputStream out;
    PrintStream p;

  	// Initiate the seed from the current time
    GregorianCalendar t = new GregorianCalendar();
    int t1 = t.get(Calendar.SECOND);
    int t2 = t.get(Calendar.MINUTE);
    int t3 = t.get(Calendar.HOUR_OF_DAY);
    int t4 = t.get(Calendar.DAY_OF_MONTH);
    int t5 = t.get(Calendar.MONTH) + 1;
    int t6 = t.get(Calendar.YEAR);
	seed = t6+70*(t5+12*(t4+31*(t3+23*(t2+59*t1))));

    if((seed%2) == 0) seed = seed-1;
    
  	initial(a);			//initializes spins 
  	
  	System.out.println("\n" + "Initial spins: ");
	//open new file	
	try{
       out = new FileOutputStream("magFlip.dat",false);
       p = new PrintStream(out);             
       p.print("\n");
       p.close();
    }
	 catch(Exception e){
           System.err.println("Error in printing to file");
     }	     

  	double accept = 0;
  	
//	System.out.println(accept + " " + i);
  	for (int i=0; i<neq; ++i){
  		iaccept = 0;
  		T = Teq;
   		flipSpin();  //exchange of spins between neighbours
  		accept = i*accept + iaccept*100.0/(lx*ly);
  		accept /= i+1;
  	}
 
    System.out.println("\n" + "After " + neq + " exchanges: ");
    try{
		out = new FileOutputStream("magFlip.dat",true);
		p = new PrintStream(out);             
        p.print("\n");
        p.println("After " + neq + " exchanges");
        p.println("Accept rate = " + accept + "%");
        double avgMag = (double)(at)/(lx*ly);	
        p.println("Average magnetization = " + avgMag);		 	              
	    p.close();
	}
	catch(Exception e){
	     System.err.println("Error in printing to file");
	}
		     
    System.out.println("Accept rate = " + accept + "%");
    
    for(int i=0; i<ntotal; ++i) {
    	iaccept = 0;
    	T = T1;
        flipSpin();
   		accept = (neq+i)*accept + iaccept*100.0/(lx*ly);
  		accept /= neq+i+1;
        if(i%nskip == 0){
      		System.out.println("\n" + "i = " + i); 
      		double avgMag = (double) (at)/(lx*ly); 
      		System.out.println("Average magnetization = " + avgMag);

			try{
		       out = new FileOutputStream("magFlip.dat",true);
		       p = new PrintStream(out);             
		       p.print("\n");
		       p.println("i = " + i);
		       p.println("Accept rate = " + accept + "%" );	
		       avgMag = (double)(at)/(lx*ly);	
        	   p.println("Average magnetization = " + avgMag);		 	              
		       p.close();
		     }
		     catch(Exception e){
		           System.err.println("Error in printing to file");
		     }

      		 System.out.println("Accept rate = " + accept + "%");
        }
    }


  	at = 0;
  	for(int i=0; i<lx; ++i)
  		for(int j=0; j<ly; ++j)
  			at += a[i][j];
  			  
  	System.out.println("Final output");

    try{
       out = new FileOutputStream("Flip.dat",true);
       p = new PrintStream(out);             
       p.print("\n");
       p.println("Final output");             
       p.close();
     }
     catch(Exception e){
           System.err.println("Error in printing to file");
     }
      
    System.out.println("Accept rate = " + accept + "%");
  }

  
//Initialize the spins
  public static void initial(int a[][]){
  	int i1, i2, j1, j2;
  	for(int i=0; i<lx; ++i)
  		for(int j=0; j<ly; ++j){
  			r1 = ranf();
  			if(r1<=fn)
   			   a[i][j] = -1;
  			else
  			   a[i][j] = 1; 
  		}
  	
  	at = 0;
  	for(int i=0; i<lx; ++i)
  		for(int j=0; j<ly; ++j)
  			at += a[i][j];
  			
 	
  	System.out.println("Excess spin = " + at );
 
 // balance spin 
  	if(at<0){
      int count = 0;
      while(count<(int)(-at/2)){
      	r1 = lx*ranf();
        int l1 = (int) r1;
        r1 = ly*ranf();
    	int l2 = (int) r1;
    	if(a[l1][l2]==-1){
    		a[l1][l2] = 1;
    		count++;	
    	}     			
      }    	
    }
    
    if(at>0){
       int count = 0;
       while(count<(int)(at/2)){
       		r1 = lx*ranf();
    	    int l1 = (int) r1;
    	    r1 = ly*ranf();
    	    int l2 = (int) r1;
    	    if(a[l1][l2]==1){
    	       a[l1][l2] = -1;
    	       count++;	
    	    }     			
    	}    	
    }
    
    at = 0;
  	for (int i=0; i<lx; ++i)
  	    for(int j=0; j<ly; ++j){
  	    	at += a[i][j];
  	    	colorSquare(i,j);	//color lattice
  	    }
  			
  	System.out.println("Excess spin = " + at);  	
 
  }
	
	
	
//////////////////////////////////////////////////////
//Flip spin
 public static void flipSpin(){
 	int i1=0, i2=0, j1=0, j2=0;
 	int nb[][] = new int[lx][ly];
 	double delH = 0;
 	double kB = 1;
 	double p; 
 
 	int daa=0, dbb=0, dab=0;
 
 	iaccept = 0;
 	//Exchange of spins 

	for(int i=lx-1; i>=0; --i)
	 	for(int j=ly-1; j>=0; --j){	
	 	    i = (int) (lx*ranf()); 	
	 	    j = (int) (ly*ranf());     
	        i1 = (lx+i-1)%lx;	//periodic boundary condition
	     	i2 = (i+1)%lx;
	     	j1 = (ly+j-1)%ly;	     
	     	j2 = (j+1)%ly;	     
	   	   
	   	  	delH = a[i1][j]+a[i2][j]+a[i][j1]+a[i][j2];
	   	 	delH *= 2.0*J*a[i][j];
	   	 	p = Math.exp(-delH/(kB*T));
		    if(p>=ranf()){
		       a[i][j] = -a[i][j];
		       iaccept += 1;
		       colorSquare(i,j);
		    }	
		}          
       
 	at = 0;
	for(int i=0; i<lx; ++i)
	    for(int j=0; j<ly; ++j)
			at += a[i][j];
 } 
 
 
 //////////////////////////////////////////////////////////////
 
 // Method to generate a uniform random number in [0,1]
// following x(i+1)=a*x(i) mod c with a=pow(7,5) and
// c=pow(2,31)-1.  Here the seed is a global variable.

  public static double ranf() {
    final int a = 16807, c = 2147483647, q = 127773,
      r = 2836;
    final double cd = c;
    int h = seed/q;
    int l = seed%q;
    int t = a*l-r*h;
    if (t > 0) seed = t;
    else seed = c+t;
    return seed/cd;
  }

 ////////////////////////////////////////////////////////////
 // given a lattice site, color that square according to the current a value:
 private static void colorSquare(int i, int j) {
    if (a[i][j] == 1) g.setColor(upColor); else g.setColor(downColor);
    g.fillRect(squareWidth*i + frameLeft, squareWidth*j + frameTop, squareWidth, squareWidth);
 }
   
}