//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//    Program ising
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

#include "cip.h"
#include "nxgraph.h"

//---- definition of fundermental constants
//---- physical constants
const double J      = 1.0;           // bonding constant
const double H      = 0.0;           // external magnetic field
const int SIZE      = 16;            // size of lattice

//---- experimental settings
// settings about temperature
const double T_INIT = 0.25;          // initial temperature
const double T_FINA = 5.00;          // final   temperature
const double T_STEP = 0.25;          // increments

// settings about sweep the lattice
const int RELAX     = 400;           // steps till relaxation
const int INTERVAL  = 10;            // steps per one measure

// settings about number of data to measure fluctuations of ...
const int SAMPLES   = 40;            // energy, magnetization
const int GROUPS    = 20;            // susceptibility, specific heat

// file names of output data
const char* fnames[4] = {
  "energy.dat", "magnet.dat", "suscep.dat", "spheat.dat"
};
// titles of each quantities
const char* titles[4] = {
  "Energy", "Magnetization", "Susceptibility", "Specific heat"
};
// file pointers
static FILE* fptrs[4];

//---- index indicators
const int ENERGY    = 0;
const int MAGNET    = 1;
const int SUSCEP    = 2;
const int SPHEAT    = 3;

//---- graphics settings
const int    WIN_WIDTH  = 256;
const int    WIN_HEIGHT = 256;

//---- declaration and definition of class Statistic
class Statistic
{
  int n;                // current total of data
  double sum1, sum2;    // current sum of data and its square
public:
  inline Statistic( void ){
    n = 0;
    sum1 = sum2 = 0.0;
  }
  //-- adds a new datum
  inline void Add( double datum ){
    n++;
    sum1 += datum;
    sum2 += datum*datum;
  }
  //-- current mean value
  inline double Mean( void ){
    return sum1/n;
  }
  //-- current variance value
  inline double Variance( void ){
    return (sum2 - sqr(sum1)/n)/(n-1);
  }
  //-- current deviation value
  inline double Deviation( void ){
    return sqrt( Variance() );
  }
  //-- current deviation of mean value, namely error value
  inline double Error( void ){
    return sqrt( Variance()/n );
  }
};


//---- declaration and difinition of class Spin
class Spin
{
  int spin;      // the value of spin, it takes +1 or -1
public:
  //-- sets this spin randomly +1 or -1
  void RandSet( void ){
    if( Urand() > 0.5 ) spin = +1; else spin = -1;
  }
  //-- represents the value of spin when this object is called
  operator int (){
    return spin;
  }
  //-- transits this spin by given probability
  void Transit( int s_sum, double prob ){
    if( Urand() < prob ) spin = +1; else spin = -1;
  }
};


//---- declaration of class Lattice
class Lattice
{
  Spin spin[SIZE][SIZE];    // spins in this lattice
  double prob[5];           // table of transition probability
public:
       Lattice( double T );
  void Sweep( int interval );
  void Measure( Statistic& Energy, Statistic& Magnet );
  void Draw( void );
private:
  inline void Neighbor( int& im, int& ip, int i );
};

//---- makes transition propability table and let spins randomly turn
Lattice::Lattice( double T )
{
  for( int s_sum=-4 ; s_sum<=+4 ; s_sum+=2 ){
    double eplus = exp( (J*s_sum+H)/T );
    prob[s_sum/2+2] = eplus/(eplus+1.0/eplus);
  }
  for( int ix=0 ; ix<SIZE ; ix++ ){
    for( int iy=0 ; iy<SIZE ; iy++ ){
      spin[ix][iy].RandSet();
    }
  }
}


//---- gets neighbors' indices of given index
inline void Lattice::Neighbor( int& im, int& ip, int i )
{
  if( i==0 ) im = SIZE-1; else im = i-1;
  if( i==SIZE-1 ) ip = 0; else ip = i+1;
}

//---- sweeps this lattice [interval] times
void Lattice::Sweep( int interval )
{
  int ix, iy, ixm, ixp, iym, iyp;
  int s_sum;
  
  while( interval-- ){
    for( ix=0 ; ix<SIZE ; ix++ ){
      Neighbor( ixm, ixp, ix );
      for( iy=0 ; iy<SIZE ; iy++ ){
        Neighbor( iym, iyp, iy );
	
        s_sum = spin[ixp][iy] + spin[ixm][iy]
              + spin[ix][iyp] + spin[ix][iym];
	
        spin[ix][iy].Transit( s_sum, prob[s_sum/2+2] );
      }
    }
  }
}

//---- evaluates energy and magnetization of this lattice
void Lattice::Measure( Statistic& Energy, Statistic& Magnet )
{
  int ix, iy, ixm, ixp, iym, iyp;
  int s_sum=0, ss_sum=0;
  
  for( ix=0 ; ix<SIZE ; ix++ ){
    Neighbor( ixm, ixp, ix );
    for( iy=0 ; iy<SIZE ; iy++ ){
      Neighbor( iym, iyp, iy );
      
      s_sum  += spin[ix][iy];
      ss_sum += spin[ix][iy] * (spin[ixm][iy] + spin[ix][iym]);
    }
  }
  
  const double energy = -J*ss_sum - H*s_sum;
  const double magnet = (double) s_sum;
  
  Energy.Add( energy );
  Magnet.Add( magnet );
}

//---- draws spins as colored square
void Lattice::Draw( void )
{
  const int dX = WIN_WIDTH/SIZE, dY = WIN_HEIGHT/SIZE;

  NXClearWindow();

  for( int ix=0 ; ix<SIZE ; ix++ ){
    for( int iy=0 ; iy<SIZE ; iy++ ){
      if( spin[ix][iy] == +1 ) NXSetColor( NX_YELLOW );
      else                     NXSetColor( NX_BLUE   );

      NXFillRectangle( ix*dX, iy*dY, dX-1, dY-1 );
    }
  }
  NXFlush();
}


//---- opens four files to save data
void OpenFiles( void )
{
  for( int i=0 ; i<4 ; i++ ){
    fptrs[i] = fopen( fnames[i], "wt" );
    if( fptrs[i] == NULL ){
      perror("Error at fopen.");
      exit(1);
    }
  }
}

//---- prints and saves each quantities in the four files
void Output( double T, Statistic quantity[] )
{
  printf("\nTemperature    = % f\n", T );

  for( int i=0 ; i<4 ; i++ ){
    printf("%-14s = % f +/- %f\n", titles[i], quantity[i].Mean(), quantity[i].Error() );
    fprintf( fptrs[i], "%f % f %f\n", T, quantity[i].Mean(), quantity[i].Error() );
  }
}

//---- closes the files
void CloseFiles( void )
{
  for( int i=0 ; i<4 ; i++ ){
    fclose( fptrs[i] );
  }
}

//---- executes experiment at given temperature
void Experiment( double T )
{
  Lattice   latt(T);
  Statistic quantity[4];

  latt.Sweep( RELAX );     // makes the lattice equilibrium

  for( int group=0 ; group<GROUPS ; group++ ){
    Statistic Energy0, Magnet0;
    
    // measures fluctuations of energy and magnetization
    for( int sample=0 ; sample<SAMPLES ; sample++ ){
      latt.Sweep( INTERVAL );
      latt.Measure( Energy0, Magnet0 );
    }
    // calculates other data per spin statistically
    quantity[ENERGY].Add( Energy0.Mean()/(SIZE*SIZE) );
    quantity[MAGNET].Add( Magnet0.Mean()/(SIZE*SIZE) );
    quantity[SUSCEP].Add( Magnet0.Variance()/T/(SIZE*SIZE) );
    quantity[SPHEAT].Add( Energy0.Variance()/sqr(T)/(SIZE*SIZE) );
    
    latt.Draw();  
  }

  Output( T, quantity ); // saves data in files
}

//---- main function
int main(void)
{
  NXOpenWindow( "Ising Model", WIN_WIDTH, WIN_HEIGHT );
  NXEnable( NX_DOUBLEBUFFER );

  OpenFiles();

  for( double T=T_INIT ; T<=T_FINA ; T+=T_STEP ){  // temperature loop
    Experiment( T );
    if( !NXNoEvents() ){
      fprintf( stderr, "Emergent stop.\n");
      break;
    }
  }

  CloseFiles();

  NXCloseWindow();
  return 0;
}