// Simulacao de Monte Carlo do modelo de trocas cineticas

// VOCE DEVE RENOMEAR O ARQUIVO ANTES DE USAR, MUDE PARA: keom.c

// INSTRUCOES PARA USUARIOS DE LINUX:

// COMPILACAO: gcc keom.c -o prog -lm

// ISSO VAI GERAR UM ARQUIVO EXECUTAVEL COM O NOME: prog 

// DEPOIS BASTA FAZER ./prog



#include <stdio.h>
#include <math.h>
#include <stdlib.h>


#define N 400  // tamanho da populacao

#define tmax 30000 // tempo maximo de simulacao (numero de passos de MC)

#define t_steady 20000 // numero de passos de MC descartados

#define cnt 10 // distancia entre medidas (NAO DEVE SER MODIFICADO)


#define nsamples 5 // numero de simulacoes independentes ("ensembles")

// gerador de numero aleatorios:
#define rfactor (long int)302875106592253 // numero "magico"
#define rmaxint (double)18446744073709551616.0 //maior inteiro p/ 64 bits


double o[N]; // opinioes

double U[N][N]; // interacoes

unsigned r0,e;

long int r;

double num,num2,num3,num4;


void main()
{
  int i,j,t,sample,cont,k,count,cont2,total;
  double mag,mag2,mag4,oi,oj,d,m,aux,mag_aux;
  double prob1,prob2,p;
  double m_ave,m2_ave,m4_ave,order_par,X,binder_cum,ratio;

  FILE *fq;

  char out_name[100];

  // abrindo o arquivo de dados:
  sprintf(out_name,"data_m_X_U_xp_N%d.dat",N); 
  fq=fopen(out_name,"w");

  r0=5454411; //semente do gerador de numeros aleatorios

  r=r0;
  for(k=0;k<100;k++) r*=rfactor;

  prob1=1.0/3.0;

  prob2=2.0/3.0;

  for(p=0.0;p<0.41;p+=0.01) //laco do parametro p
    {                        //p: fracao de interacoes negativas

      m_ave=0.0;
      m2_ave=0.0;
      m4_ave=0.0;
      
      order_par=0.0;
      X=0.0;
      binder_cum=0.0;



      for(sample=0;sample<nsamples;sample++) // laco de simulacoes diferentes
	{

	  mag=0.0;
	  mag2=0.0;
	  mag4=0.0;
	  total=0;
	  count=0;
	   
	  // condicao inicial:
         
	  for(i=0;i<N;i++)
	    {
	      // opinioes iniciais:
	      // comecando com 1/3 de cada opiniao
	      r*=rfactor;
	      num=((unsigned long)r)/rmaxint;
	      
	      if(num<prob1) o[i]=1.0;
	      else if(num<prob2) o[i]=-1.0;
	      else o[i]=0.0;
	      
	      // interacoes entre todos os pares de agentes:

	      for(j=0;j<N;j++)
		{
		  if(j==i) U[i][j]=0.0;
		  else
		    {
		      r*=rfactor;
		      num2=((unsigned long)r)/rmaxint;

		      // U=-1 com probabilidade p
		      // U=+1 com probabilidade 1-p
		      
		      if(num2<p) U[i][j]=-1.0;
		      else U[i][j]=1.0;
		
		    }
		}

	      
	    }

  
	  //############################################################
	  
	  //                        Dinamica

	  //############################################################
	   

	  for(t=0;t<tmax;t++) //laco de tempo
	    {	
	      
	      cont2=0;
	      
	      do
		{
		  
		  
		  
		  // Escolhemos aleatoriamente 2 agentes: "i" e "j"
		  r*=rfactor;
		  num3=((unsigned long)r)/rmaxint;
		  
		  i=((int) N*num3);
		  
		  do
		    {
		      r*=rfactor;
		      num4=((unsigned long)r)/rmaxint;
		      
		      j=((int) N*num4);
		    }while(j==i);
		  
		  
		  // Interacao:
		  
		  oi=1.0*o[i]+1.0*U[i][j]*o[j];

		  // garantindo que nao tenhamos o>1 ou o<-1:
		  if(oi>1.0) oi=1.0;
		  if(oi<-1.0) oi=-1.0;
		  
		  o[i]=oi;
		  
		  
		  cont2++;
		  
		}while(cont2<N); //1 passo de tempo corresponde a N interacoes


	      if(t>t_steady) //O sistema atingiu o estado estacionario?
		{
		  count++;
		  if(count==cnt)
		    {
		      count=0;
		      
		      m=0.0;

		      /// calculo do parametro de ordem:
		      for(k=0;k<N;k++) m+=o[k];
		      
		      mag_aux=fabs(m);
		      
		      mag+=mag_aux;
		      mag2+=( mag_aux*mag_aux );
		      mag4+=( mag_aux*mag_aux*mag_aux*mag_aux );
		      
		      total++;  //conta o numero de pontos que serao medidos
		      
		      
		    }

		}

	    } // fim do laco do tempo

	  
	  // medias temporais
	  mag/=total;  
	  mag2/=total;
	  mag4/=total;
	  
	  
	  m_ave+=mag;  
	  m2_ave+=mag2;
	  m4_ave+=mag4;
	  
	} // fim do laco de simulacoes diferentes
      

      // medias sobre simulacoes diferentes
      m_ave/=nsamples;  
      m2_ave/=nsamples;
      m4_ave/=nsamples;
      
      order_par=((double) 1.0*m_ave/N);
      
      ratio=((double) m4_ave/(3.0*m2_ave*m2_ave));
      
      binder_cum=((double) 1.0-1.0*ratio);
      
      X=( m2_ave-(m_ave*m_ave) )/(N);

      printf("m(p=%3.2lf)=%12.10lf\n",p,order_par);

      // imprime os dados no arquivo
      // 1a coluna: p
      // 2a coluna: parametro de ordem
      // 3a coluna: "susceptibilidade"
      // 4a coluna: cumulante de Binder
      
      fprintf(fq,"%3.2lf\t%12.10lf\t%12.10lf\t%12.10lf\n",p,order_par,X,binder_cum);
      
      fflush (fq);
      
      
    } // fim do laco para o parametro "p"
  
  
  fclose(fq);
  
}
    
