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

#define MAXN 1000000

double x[MAXN];       // Positions of particles are globally defined 
double deltaT = 0.1;  // Time step
int Ntime = 100;      // Number of time steps

int init(int);
double averageDistance(int); 
int update(int);
double move(double, double);

int main(int argc, char *argv[]) {
  int numP, i;
  double avgeDist, data[2];

  if(argc != 2) {
    printf("Usage: %s  number_of_particles\n", argv[0]);
    exit(EXIT_FAILURE);
  }
  numP = atoi(argv[1]);
 
  /* Initialize positions */
  init(numP);

  /* Average position after initialization */
  avgeDist = averageDistance(numP);
  data[0] = avgeDist;

  /* Loop over timesteps for update using all particles */
  for(i=0; i<Ntime; i++){
    update(numP);
  }

  /* Average position after the update() */
  avgeDist = averageDistance(numP);
  data[1] = avgeDist;

  printf("Average distance between particles at timestep 0 is %f\n", data[0]);
  printf("Average distance between particles at timestep %d is %f\n", Ntime, data[1]);

  return 0;
}

/*******************************************************************/
/* Initialize the positions of particles to be evenly distributed. */
/*******************************************************************/
int init(int n) {
  int i;
  for(i=0; i<n; i++){
    x[i] = (double)i / (double) n;
  }
  return 0;
}

/****************************************************************************/
/* Update the positions of particles using forces from all other particles. */
/****************************************************************************/
int update(int n) {
  int i,j;
  double totalForce[MAXN];

  for(i=0; i<n; i++) {
    totalForce[i] = 0.0;
    for(j=0; j<n; j++) {
      if (j < i) totalForce[i] -= 1.0 / pow(x[i]-x[j], 2);
      if (j > i) totalForce[i] += 1.0 / pow(x[i]-x[j], 2);
    }
  }

  for(i=0; i<n; i++) {
    x[i] = move(x[i], totalForce[i]);
  }
  return 0;
}

/*************************************************************/
/* Find the average distance between two adjacent particles. */
/*************************************************************/
double averageDistance(int n) {
  int i;
  double distance;

  /* Loop over all particles */
  distance = 0.0;
  for(i=1; i<n; i++) {
    distance += x[i] - x[i-1];
  }
  /* Compute average distance */
  distance = distance / (n-1);
  return distance;
}

/********************************************************/
/* Update the position of a particle using a simplified */
/* Newton's law of motion (ignore mass and velocity).   */
/********************************************************/
double move(double xinit, double force) {
  double xfinal;
  xfinal = xinit + 0.5 * force * deltaT * deltaT;
  return xfinal;
}