/* Vicsek model simulation -- University of Warwick, 2013 */

/* compile using: gcc vicsek.c -o vicsek */
/* run using: ./vicsek [size] [N] [range] [noise] [timesteps] [stepskip] [seed] */

/* BEGIN HEADER */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

/* global variables */
unsigned int STATE[32];
unsigned long int seed;

/* prototypes */
void parameters(int nump, char *param[], double *size, int *N, double *range, double *noise, int *timesteps, int *stepskip);
/* END HEADER */

//=============================
int main(int argc, char *argv[])
{
  char outputfile[80],flockname[80],ordername[80];
  int N, timesteps, stepskip, neighbours, count;
  double size, range, noise, dx, dy, dt, dist, mx, my, order, av_order;
  double *posx, *posy, *theta, *D_theta;
  register int i, j, k;
  FILE *flock_output, *order_output;

  /* Check parameters */
  parameters(argc, argv, &size, &N, &range, &noise, &timesteps, &stepskip);
  dt = 0.25; // timestep
  av_order = 0.0; count = 0; // variables for calculating the average order

  sprintf(flockname,"flock_%d.dat",(int) N); // add flock size to file name
  flock_output = fopen(flockname,"a"); // open file
  sprintf(ordername,"order_%d.dat",(int) N);
  order_output = fopen(ordername,"a");
 
  posx = malloc(N*sizeof(int*)); // allocate memory
  posy = malloc(N*sizeof(int*));
  theta = malloc(N*sizeof(int*));
  D_theta = malloc(N*sizeof(int*));
  
  // initialise a random configuration for the flock
  for (i=0; i<N; i++) { 
    posx[i] = size*drand48();
    posy[i] = size*drand48();
    theta[i] = 2.0*M_PI*drand48();
  }
  
  /* Simulate */
  for (i=1; i<=timesteps; i++) {  // for each timestep
    
    order = 0.0; mx = 0.0; my = 0.0;
    
    // run over the flock
    for (j=0; j<N; j++) {
      // propagate
      posx[j] += dt*cos(theta[j]); 
      if (posx[j] > size) {posx[j] -= size;} 
      if (posx[j] < 0.0) {posx[j] += size;}
      posy[j] += dt*sin(theta[j]); 
      if (posy[j] > size) {posy[j] -= size;}
      if (posy[j] < 0.0) {posy[j] += size;}
      // order parameter
      mx += cos(theta[j]); my += sin(theta[j]);
      // align -- expensive O(N^2) algorithm
      neighbours = 0; D_theta[j] = 0.0;
      for (k=0; k<N; k++) { 
	dx = posx[k]-posx[j]; dy = posy[k]-posy[j];
	if (dx > size/2.0) {dx -= size;} 
	if (dx < -size/2.0) {dx += size;} 
	if (dy > size/2.0) {dy -= size;} 
	if (dy < -size/2.0) {dy += size;} 
	dist = sqrt(dx*dx+dy*dy);
	if (dist < range) {
	  neighbours++;
	  D_theta[j] += theta[k]-theta[j];
	}
      }
      if (neighbours != 0) D_theta[j]/=((double) neighbours);
      order = sqrt(mx*mx+my*my)/((double) N);
    }

    // update
    for (j=0; j<N; j++) {
      theta[j] += D_theta[j] + noise*M_PI*(drand48()-1.0); 
      if (theta[j] > 2.0*M_PI) {theta[j] -= 2.0*M_PI;}
      if (theta[j] < 0.0) {theta[j] += 2.0*M_PI;}
    }
    
    // output flock
    if (i % stepskip == 0) {
      for (j=0; j<N; j++) {
      	fprintf(flock_output,"%.10f\t%.10f\t%.10f\n",posx[j],posy[j],theta[j]);
      }
      fprintf(order_output,"%d\t%.10f\n",i,order);
      count++;
      av_order = av_order*(((double) count-1)/((double) count)) + order/((double) count);
    }
  } // end sweep

  printf("the average degree of alignment is %.6f\n",av_order);
  
  free(posx); // free memory
  free(posy); // free memory
  free(theta); // free memory
  free(D_theta); // free memory
  
  fclose(flock_output); // close file
  fclose(order_output); // close file
  
  return EXIT_SUCCESS;
} // end main
//==============================

/* Check user-provided parameters */
void parameters(int nump, char *param[], double *size, int *N, double *range, double *noise, int *timesteps, int *stepskip)
{	

  if (nump<8) {
    printf("\nUsage: %s [size] [N] [range] [noise] [timesteps] [stepskip] [seed]\n\n",param[0]);
    printf("\t[size]      is the dimension of each side of the periodic box.\n");
    printf("\t[N]         is the size of the flock.\n");
    printf("\t[range]     is the interaction range.\n");
    printf("\t[noise]     is the strength of the noise.\n");
    printf("\t[timesteps] is the number of timesteps to run.\n");
    printf("\t[stepskip]  is the frequency to output data.\n");
    printf("\t[seed]      is the random number generator seed (positive integer).\n");
    exit(EXIT_FAILURE);
  }
  
  *size = atoi(param[1]);
  if (*size<2.0) {
    printf("\nThe size of the model must be at least 2x2.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *N = atof(param[2]);
  if (*N<1) {
    printf("\nThe number of birds must be at least 1.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *range = atof(param[3]);
  if (*range<0.0 || *range>*size) {
    printf("\nThe range of interaction must be between 0.0 and the periodic size.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *noise = atof(param[4]);
  if (*noise<0.0) {
    printf("\nThe noise strength must be positive.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *timesteps = atoi(param[5]);
  if (*timesteps<1) {
    printf("\nYou must run at least 1 timestep.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *stepskip = atoi(param[6]);
  if (*stepskip<1 || *stepskip>*timesteps) {
    printf("\nYou must output data.\n\n");
    exit(EXIT_FAILURE);
  }
  
  seed=strtoul(param[7],NULL,10);
  if (seed<=0) {
    printf("\nThe random number seed must be a positive integer.\n\n");
    exit(EXIT_FAILURE);
  }
  srand48(seed);

  return;
} // end parameters