//*************Ising Model Viewing***********************//

/* ISING MODEL
 *
 * Michael Stefferson
 * Independent Study PHYS 492
 * Advisor: D. Toussaint
 * January 21, 2011
 *
 * Compiling: gcc -o isingmodel_view isingmodel_view.c /usr4/mws2/ISING/animation/libanimation.a -lX11 -lm
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include "/usr4/mws2/ISING/animation/animation.h"

#define MAX_PARTICLES 64

double probab(double E);
/*
int n_particles;

float x[MAX_PARTICLES], y[MAX_PARTICLES], spin[MAX_PARTICLES];
float LX, LY;

float dt, mytime, tmax;

int n_particles;
*/

int main(void)
{

int n_particles;

float x[MAX_PARTICLES], y[MAX_PARTICLES], spin[MAX_PARTICLES];
float LX, LY;

float dt, mytime, tmax;


  int atoms1D = 8, elements;

  // int atoms1D = sqrt(atoms1D);
  int particle[atoms1D][atoms1D], part_position[atoms1D][atoms1D];
  int i,j, k, k_end, R, l, t, b;

  //If I want to vary individual parameters
  double kb, J, T, B;
  int spinup = 1, spindown = -1, tot_spinup, tot_spindown, mag_spins, magnet_count;
  double probE,energy, probup, R, aver_magnet, aver_energy;
  int iter, counter, iter_end, iter_count;
  int seed = 4;
  int method;

  //animation
  int anim_i;



  srand48(seed);

  //fprintf(stderr,"What do you want from me?\n1: matrices\n2: interations vs. max spin percent\n3: M vs. B\n4: M vs T (~1/B)\n5: Average E vs  B\n6: Average E vs. T\n");

  // scanf("%d", &method);


 tot_spinup = 0;
 tot_spindown = 0;
 elements = 0;
 mag_spins = 0;


 n_particles = 64;
 LX =  1.0;
 LY = 1.0;
 dt = 0.005;

 x[0] = 0; y[0] = 0;

 //position stuff
 for(i = 1; i < n_particles; i++){
     if( i % atoms1D == 0 ){
       i == 0;
     }
     x[i] = (double)i * LX / atoms1D;
   }
 for(i = 1; i < n_particles; i++){
   if( i % atoms1D == 0 ){
     i == 0;
   }
   y[i] = (double)i *  LY / atoms1D;
 }


 /*
 //position stuff
 for(i = 1; i < n_particles; i++){
   for(j = 1; j < n_particles; j++){
posi
 */


//animation stuff
//animation_start( 0.0, LX, 0.0, LY );
//anim_i = add_many_particles(n_particles, x, y);
//set_clock_tick(dt * 10.0);
tmax = 500;


  //initialize spins
  for(i = 0; i < atoms1D; i++){
    for(j = 0; j < atoms1D; j++){
      tot_spinup ++;
      elements ++;
      particle[i][j] = spinup;
      mag_spins += particle[i][j];
      // printf("mag spins: %d", mag_spins);
    }
  }



  counter = 0 ;

  for(mytime = 0; mytime <= tmax; mytime += dt){
  //right, left, above, below. periodic B.C.

    tot_spinup = 0;
    tot_spindown = 0;
    elements = 0;
    mag_spins = 0;


    for(i = 0; i < atoms1D; i++){
      for(j = 0; j < atoms1D; j++){
    //energy/kbT for the ith, jth component. Average the neighbors.

    R = i + 1;
    l = i - 1;
    t = j + 1;
    b = j - 1;

      //left edge
    if( i == 0 )
    {
      l = atoms1D - 1;
    }
    //bottom edge
    if( j == 0 )
      {
        b = atoms1D - 1;
      }

    //right
    if( i == (atoms1D - 1))
      {
        R = 0;
      }

    //top
    if( j == (atoms1D - 1) )
      {
        t = 0;
      }


    probE = -B  *( (double)particle[R][j] + (double)particle[l][j] + (double)particle[i][t] +  (double)particle[i][b] );

    probup = probab(probE);

      //  printf("E: %lf probup: %lf\n",E, probup);

      //monte carlo
      R = drand48();

      //  printf("E: %lf probup: %lf R: %lf \n",E, probup, R);

      if(R <= probup){
    particle[i][j] = spinup;
      }
      if(R > probup){
    particle[i][j] = spindown;
      }
    mag_spins += particle[i][j];
    //  printf("particle: %d mag_spins: %d\n", particle[i][j], mag_spins);
      }
    }// end spin determine LOOP

    i = 0; j = 0;
    for(k = 0; k < n_particles; k++){
      printf("%d %d\n", i,j);
      spin[k] = particle[i][j];
      printf("k: %d  spin:%d particle[%d][%d]: %d %d %d\n", k, spin[k], i, j, particle[i][j], i, j);
      i ++;
      if(i % atoms1D == 0){
        j ++;
        i = 0;
      }
    }

    // move_many_particles(n_particles, anim_i, x, y);


    //animation

    for(i = 0; i < atoms1D; i++){
      if(spin[k] == spinup){
    printf("c %f %f 1\n", x[i], y[i]);
      }

  if(spin[k] == spindown){
    printf("d %f %f 1\n", x[i], y[i]);
      }
    }
    printf("F/n");
  }//end time LOOP

  return 0;
}//end main