/codigo

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

// programa particulas boids para CC periodica. 

#define N 500 //numero de particulas
#define NUM 4 // numero de predadores
#define T 5000 //numero de passos
#define C 0.05 //constante de coesao
#define S 0.3 // constante de separaçao
#define M 0.15 // media das velocidades
#define B 0.6 // constante de repulsao da barreira
#define ini 0.1 // valor maximo das componentes da velocidade
              // valor maximo da velocidade (sqrt(2)*ini)
#define E 10 // valor de excesso que a barreira trava os BOIDS
#define EE 0.1 // valor constante barreira nas fronteiras
#define PI 3.1415926538979 // aprox de pi
#define I 2 // Quantidade de passos passados para plotar no gif
 
double norm(double x, double y){ //funçao norma
  return sqrt(x*x + y*y);
}


main(){
  
  double r[N][2], v[N][2], epsilon=0.01, R=20, rmin=3, erre=12, ERRE=10, dt = 0.1;
  double random;
  double x=0., y=0., L=100, vmx=0, vmy=0, vmed=0, vmedl =0;
  double bx[NUM], by[NUM], vx[NUM], vy[NUM]; //barreira
  int i=0, j=0, k=0, l=0, m=0, t=0, a=0, b=0;
  int cont=0, kont=0; //conta o numero de vizinhos proximos da particula
  int vizinhos[N][2]; //nao deixa a velocidade explodir nem tender a 0


  //**********Inicializar Particulas (posicao e velocidade)**********//
  
  printf("set terminal gif animate delay 1.0E-5\n");// Produz arquivo .gif
  printf("set output 'boidsbar.gif'\n");            //
  printf("set size square\n");                            // caso rode com
  printf("set title 'Tempo t = 0.0'\n");                  // |gnuplot apa-
  printf("set xrange [0:%lf]\n", L);                      //  recem os graficos.
  printf("set yrange [0:%lf]\n", L);                      // 
  printf("plot '-' w p pt 3 title sprintf(\"BOIDS\")\n"); //
  
  for(i=0;i<N;i++){
    r[i][0] = L*rand()/RAND_MAX;            //
    r[i][1] = L*rand()/RAND_MAX;            // Inicializa as pos e vel
    v[i][0] = ini*(2.*rand()/RAND_MAX - 1); // das particulas aleatoriamente
    v[i][1] = ini*(2.*rand()/RAND_MAX - 1); //
    printf("%lf %lf \n", r[i][0], r[i][1]); // 
  }

  printf("e \n"); // Final da impressao dos BOIDS
  
  
  //**********Dinamica das Particulas**********//

  
  a=0;
  do{  //inicio do passo de tempo
    t++;  //passo de tempo
    if(t%I == 0){
      printf("set size square\n");                    //
      printf("set title 'Tempo t = %.1lf'\n", t*dt);  // 
      printf("set xrange [-0:%lf]\n", L);       // Caso rode com |gnuplot 
      printf("set yrange [-0:%lf]\n", L);       // sao gerados os graficos
      printf("plot '-' w p pt 3 title sprintf(\"BOIDS\"), '-' w p pt  7 title sprintf(\"Predadores\")\n"); // 
      }
    
      vmed = 0;  // velocidade media das particulas
      
    
    //**********Inicializar Barreiras (Predadores)**********//
    
    if(0 == t%700 && a < NUM){ // Predadores Adicionados de 500 em 500 passos
      if(a%4 == 0){ // Ou os predadores estao posicionados em cima
	bx[a] = L*rand()/RAND_MAX; // Inicializar barreira na componente x
	by[a] = L;                 // Inicializar barreira na componente y
	vx[a] = -1.3*ini*(rand()/RAND_MAX+0.5); //Inicializar velocidade em x
	vy[a] = -1.3*ini*(rand()/RAND_MAX+0.5); //Inicializar velocidade em y
      }else if(a%4 == 1){ // Ou predadores sao adicionados na lateral
	bx[a] = 0;                 // Inicializar barreira na componente x
	by[a] = L*rand()/RAND_MAX; // Inicializar barreira na componente y
	vx[a] = 1.3*(ini*rand()/RAND_MAX+0.5);   //Inicializar velocidade em x
	vy[a] = -1.3*ini*(rand()/RAND_MAX+0.5); //Inicializar velocidade em y
      }else if(a%4 == 2){
	bx[a] = L*rand()/RAND_MAX; // Inicializar barreira na componente x
	by[a] = 0; // Inicializar barreira na componente y
	vx[a] = 1.3*ini*(rand()/RAND_MAX+0.5);   //Inicializar velocidade em x
	vy[a] = 1.3*ini*(rand()/RAND_MAX+0.5); //Inicializar velocidade em y
      }else if(a%4 == 3){
      	bx[a] = L; // Inicializar barreira na componente x
	by[a] = L*rand()/RAND_MAX; // Inicializar barreira na componente y
	vx[a] = -1.3*ini*(rand()/RAND_MAX+0.5);//Inicializar velocidade em x
	vy[a] = 1.3*ini*(rand()/RAND_MAX+0.5); //Inicializar velocidade em y
      }
      a++;  // Vai para o proximo contador
    }
    
    for(j=0;j<N;j++){  // dinamica de cada particula
      
      cont=0; // contador para coesao das particulass
      x = 0; // media das particulas vizinhas para coesao
      y = 0; // media das particulas vizinhas para coesao
      
      kont = 0; // contador para media das velocidades
      vmx = 0; // media das veloc. das particulas
      vmy = 0; // media das veloc. das particulas
      
      
      //**********Alinhamento**********//
      
      
      for(k=0;k<N;k++){ // laço para todas as particulas
	if(k == j){        // se eh a mesma particula nao faz nada
	} else if(norm(r[j][1]-r[k][1],r[j][0]-r[k][0])<erre){ 
	  kont++;              // se eh outra particula acrescenta no contador
	  vmx = vmx + v[k][0]; // acrescenta a velocidade na componente x
	  vmy = vmy + v[k][1]; // acrescenta a velocidade na componente y
	}
      }
      
      if(kont > 1){
	v[j][0] = (1-M)*v[j][0] + M*vmx/kont; // media para a velocidade
	v[j][1] = (1-M)*v[j][1] + M*vmy/kont; // dos vizinhos e a propria
      }
      
      
      //**********Coesao Entre Particulas e CM**********//
      
      
      for(m=0;m<N;m++){ // laço das particulas 
	if(norm(r[j][1]-r[m][1],r[j][0]-r[m][0])<R){// se as particulas tiverem
	                                            // a uma distancia menor q
	                                            // "R"
	  cont++;          // acrescenta no contador 
	  x = x + r[m][0]; // acrescenta na posicao das particulas na comp. x
	  y = y + r[m][1]; // acrescenta na posicao das particulas na comp. y
	}	 
      }
      
      if(cont>1){ // se contador for cont > 1
	x = x/cont;     // media aritmetica na componente x       
	y = y/cont;     // media aritmetica a componente y
	
	v[j][0] = v[j][0] - C*(r[j][0] - x)/L; // atualiza velocidade c poten-
	                                       // cial de mola na componente x
	v[j][1] = v[j][1] - C*(r[j][1] - y)/L; // atualiza velocidade c poten- 
      }                                        // cial de mola na componente y
      
      
      //**********Separacao Entre Particulas**********//
      
      for(l=0;l<N;l++){ //laço para particulas
	if(l==j){  // se a particula for ela mesma nao faz nada 
	  
	}else{     // se nao é ela entao
	  
	  if(norm(r[j][0]-r[l][0],r[j][1]-r[l][1])<rmin){ // se as particulas
	                              // estao a uma distancia menor que  "rmin"
	    v[j][0] = v[j][0] + S*(r[j][0] - r[l][0])/pow(norm(r[l][0]-r[j][0],r[l][1]-r[j][1]),2); // Atualiza velocidade c potencial na componente x
	    v[j][1] = v[j][1] + S*(r[j][1] - r[l][1])/pow(norm(r[l][0]-r[j][0],r[l][1]-r[j][1]),2); // Atualiza velocidade c potencial na componente y
	  }
	}
      }
      
      
      //**********Velocidade Pela Densidade de Vizinhos**********//
      
      if(cont  + kont < 0.1*N){   // se o numero de vizinhos for menor que 0.1*N
	vizinhos[j][0] = 1;       // chave 1 fecha
	vizinhos[j][1] = 0;       // chave 2 abre
      }else{
	vizinhos[j][0] = 0;       // chave 1 abre
      }
      
      if(norm(v[j][0],v[j][1]) > 0.5*ini){   // se a velocidade media do passo anterior for maior
	                                    // que 0.7*ini
	vizinhos[j][1] = 0; // chave 2 "abre"
      }
      
      if(vizinhos[j][0] == 0){        // se chave  1 aberta 
	if(vizinhos[j][1] == 1){      // se chave 2 fechada  entao nao faz nada
	  
	}else if(cont +kont >0.5*N){ // se vizinhos > 0.5*N entao 
	  v[j][0] = 0.8*v[j][0];     // velocidade reduz 20% 
	  v[j][1] = 0.8*v[j][1];     // velocidade reduz 20%
	  
	}else if(cont + kont > 0.4*N){ // se vizinhos 0.5> e >0.4 entao
	  v[j][0] = 0.85*v[j][0];      // velocidade reduz 15%
	  v[j][1] = 0.85*v[j][1];      // velocidade reduz 15%
	  
	}else if(cont + kont > 0.3*N){ // se vizinhos 0.4> e >0.3 entao
	  v[j][0] = 0.9*v[j][0];       // velocidade reduz 10%
	  v[j][1] = 0.9*v[j][1];       // velocidade reduz 10%
	  
	}else if(cont + kont > 0.2*N){ // se vizinhos 0.3> e >0.2 entao
	  v[j][0] = 0.95*v[j][0];      // velocidade reduz 5%
	  v[j][1] = 0.95*v[j][1];      // velocidade reduz 5%
	}
	vizinhos[j][1] = 1; // chave 2 fecha
      }
      
      
      //**********Ruido Particulas**********//
      if(rand()%2 == 0){ // criterio para decidir a fase
	random = - rand()/RAND_MAX; // fase
      }else{
	random = rand()/RAND_MAX; // fase
      }
      v[j][0] = cos(PI*random)*v[j][0] - sin(PI*random)*v[j][1];  // rotaçao
      v[j][1] = sin(PI*random)*v[j][0] + cos(PI*random)*v[j][1]; // rotaçao
      
      
      
      //Barreira puntual no centro
      
      for(b=0;b<NUM;b++){
	if(norm(r[j][0] -bx[b], r[j][1] - by[b]) < ERRE){//se a particula ta a
	                                                 //menos de "ERRE" entao
	  
	  if(r[j][0] < bx[b] && r[j][1] < by[b]){ // se estiver atras e abaixo
	    v[j][0] = v[j][0] - B*norm(r[j][0] - bx[b], r[j][0] - bx[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	    v[j][1] = v[j][1] - B*norm(r[j][1] - by[b], r[j][1] - by[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	  
	  }else if(r[j][0] > bx[b] && r[j][1] < by[b]){ // se estiver a frente e abaixo
	    v[j][0] = v[j][0] + B*norm(r[j][0] - bx[b], r[j][0] - bx[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	    v[j][1] = v[j][1] - B*norm(r[j][1] - by[b], r[j][1] - by[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	  
	  }else if(r[j][0] < bx[b] && r[j][1] > by[b]){ // se estiver abaixo e acima
	    v[j][0] = v[j][0] - B*norm(r[j][0] - bx[b], r[j][0] - bx[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5);// atualiza velocidade com potencial analogo "gravitac"
	    v[j][1] = v[j][1] + B*norm(r[j][1] - by[b], r[j][1] - by[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5);// atualiza velocidade com potencial analogo "gravitac"
	  
	  }else if(r[j][0] > bx[b] && r[j][1] > by[b]){// se estiver a frente e acima
	    v[j][0] = v[j][0] + B*norm(r[j][0] - bx[b], r[j][0] - bx[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	    v[j][1] = v[j][1] + B*norm(r[j][1] - by[b], r[j][1] - by[b])/pow(norm(r[j][0] - bx[b], r[j][1]- by[b]),1.5); // atualiza velocidade com potencial analogo "gravitac"
	  }
	}
      }
    
      //*********Um Pouco de Barreira nas Fronteiras**********//
      
      if(r[j][0] - L < E){
	v[j][0] = v[j][0] - EE*exp(-pow((L-r[j][0])/E,2));// Atualiza velocidade c 
	// potencial gaussiano na componente x
      }
      if(r[j][0] < E){
	v[j][0] = v[j][0] + EE*exp(-pow(r[j][0]/E,2)); // Atualiza velocidade c pot
	//  gaussiano na componente x
      }
      if(r[j][1] - L < E){
	v[j][1] = v[j][1] - EE*exp(-pow((L-r[j][1])/E,2));// Atualiza velocidade c pot
	//  gaussiano na componente y
      }
      if(r[j][1] < E){
	v[j][1] = v[j][1] + EE*exp(-pow(r[j][1]/E,2)); // Atualiza velocidade c pot
	//  gaussiano na componente y
      }
      
      //**********Passo das particulas com CC periodica**********//
      
      if(L-r[j][0] < epsilon){ // final do dominio na componente x 
	r[j][0] = L - r[j][0]; // recomeça no inicio do dominio em x 
      }else if(r[j][0] < 0){    // inicio do dominio na coomponente x
	r[j][0] =  L + r[j][0]; // recomeça no finaldo dominio em x
      }
      
      r[j][0] = r[j][0] + v[j][0]*dt; // atualiza a posiçao da particula em x
      
      if(L-r[j][1] < epsilon){  // final do dominio em y
	r[j][1] =  L - r[j][1]; // recomeça no inicio do dominio em y
      }else if(r[j][1] < 0){   // inicio do dominio em y
	r[j][1] = L + r[j][1]; // recomeça no final do dominio em y
      }
      
      r[j][1] = r[j][1] + v[j][1]*dt; // atualiza a posiçao da particula em y

      if(t%I == 0){
	printf("%lf %lf \n", r[j][0], r[j][1]); //imprime a posiçao da particula
      }
      vmed = vmed + norm(v[j][0],v[j][1]); // contador para velocidade media
    }

    if(t%I == 0){
      printf("e\n"); // final do grafico para |gnuplot 
    }
    
    //*********Passo da Barreira (Predador)*********//
    
    
    for(b=0;b<a;b++){
      
      if(L - bx[b] < epsilon){ // final do dominio em x
	bx[b] = L - bx[b];            // recomeça no inicio do dominio em x
      }else if(bx[b] < epsilon){ //inicio do dominio em x
	bx[b] =  L + bx[b];         // recomeça no final do dominio em x
      } 
      
      bx[b] = bx[b] + vx[b]*dt; // atualiza posiçao da barreira em x
      
      
      if(L - by[b] < epsilon){ // final do dominio em y 
	by[b] =  L - by[b];           // recommeça no inicio do dominio em y
      }else if(by[b] < epsilon){  // inicio do dominio em y
	by[b] = L + by[b];           // recomeça no final do dominio em y
      }
      
      by[b] = by[b] + vy[b]*dt; // atualiza posiçao da barreira em y

      if(t%I == 0){
	printf("%lf %lf\n",bx[b], by[b]); // imprime a barreira
      }
    }
     if(t%I == 0){
       printf("e \n"); // final do grafico no tempo t
       printf("\n\n");
     }
    vmedl = vmed/N; // media da velocidade das particulas
  }while(t<T);
  
  printf("\n\n");
}