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

#define MAX_ELK 200
FILE *fp;

void main()
{

  double x[MAX_ELK][2],v[MAX_ELK][2],f[MAX_ELK][2];
  int    n,nsteps,seed,noutput,nalive,step,i;
  double radius,vsum,distsq,radsq;
  void   initialize_elk(),force_elk(),move_elk(),dump();

/* User inputs. */

  n = MAX_ELK + 1;
  while (n > MAX_ELK) {
    printf("Enter number of elk to simulate: \n");
    scanf("%d", &n);
  }

  printf("Enter radius of wilderness: \n");
  scanf("%lf", &radius);

  printf("Enter number of time steps: \n");
  scanf("%d", &nsteps);

  printf("Enter random number seed: \n");
  scanf("%d", &seed);

  noutput = n + 1;
  while (noutput > n) {
    printf("Enter which elk to print out: \n");
    scanf("%d", &noutput);
  }

/* Initialize the elk positions and velocities. */

  initialize_elk(n,x,v,radius,seed);

/* Dump initial snapshot to file. */

  dump(n,x,0,radius);

/* Timestep loop:
    Compute force on each elk.
    Move elk under influence of forces.
    Dump snapshot of coords.
*/
  
  for (step = 1; step <= nsteps; step++) {
    force_elk(n,x,f,radius);
    move_elk(n,x,v,f,radius);
    dump(n,x,step,radius);
  }

/* Close snapshot file. */

  dump(n,x,-1,radius);

/* Compute output. */

  radsq = radius * radius;
  nalive = 0;
  vsum = 0.0;

  for (i = 0; i < n; i++) {
    distsq = x[i][0]*x[i][0] + x[i][1]*x[i][1];
    if (distsq <= radsq) {
      nalive++;
      vsum += sqrt(v[i][0]*v[i][0] + v[i][1]*v[i][1]);
    }
  }

  if (nalive > 0) vsum = vsum/nalive;

/* Print output. */

  printf("# of Elk alive = %d\n",nalive);
  printf("Average speed of live elk = %g\n",vsum);
  printf("Coords of elk %d = %g %g\n",noutput,
	 x[noutput-1][0],x[noutput-1][1]);

}


/* Initialize the position and velocity of each elk. */

void initialize_elk(n,x,v,radius,seed)
  int    n;
  double x[MAX_ELK][2],v[MAX_ELK][2];
  double radius;
  int    seed;
{

  int    i,iflag;
  double radsq,distsq,sigma;
  double rand1();

  radsq = radius * radius;
  sigma = sqrt(3.1415926/n) * radius;

/* generate random coordinates until elk is inside Gila
   generate random velocity */

  for (i = 0; i < n; i++) {
    iflag = 1;
    while (iflag) {
      x[i][0] = (2.0 * rand1(&seed) - 1.0) * radius;
      x[i][1] = (2.0 * rand1(&seed) - 1.0) * radius;
      distsq = x[i][0]*x[i][0] + x[i][1]*x[i][1];
      if (distsq <= radsq) iflag = 0;
    }
    v[i][0] = (rand1(&seed)-0.5) * sigma;
    v[i][1] = (rand1(&seed)-0.5) * sigma;
  }
}


/* Compute the forces acting on each elk. */

void force_elk(n,x,f,radius)
  int    n;
  double x[MAX_ELK][2],f[MAX_ELK][2];
  double radius;
{

  int    i,j;
  double epsilon,sigma,sigmasq,radsq,distsq;
  double delx,dely,sig2r2,force;

  epsilon = 1.0;
  sigma = sqrt(3.1415926/n) * radius;
  sigmasq = sigma * sigma;
  radsq = radius * radius;

/* Zero out forces. */

  for (i = 0; i < n; i++) {
    f[i][0] = 0.0;
    f[i][1] = 0.0;
  }

/* N^2 double loop to accumulate forces.
   Elk outside of Gila do not contribute. */

  for (i = 0; i < n; i++) {
    distsq = x[i][0]*x[i][0] + x[i][1]*x[i][1];
    if (distsq <= radsq) {
      for (j = 0; j < n; j++) {
	distsq = x[j][0]*x[j][0] + x[j][1]*x[j][1];
	if (distsq <= radsq && i != j) {
	  delx = x[i][0] - x[j][0];
	  dely = x[i][1] - x[j][1];
	  distsq = delx*delx + dely*dely;
	  if (distsq <= sigmasq) {
	    force = 2.0*epsilon/sigmasq;
	  }
	  else {
	    sig2r2 = sigmasq / distsq;
	    force = 4.0*epsilon/distsq * sig2r2 * (sig2r2-0.5);
	  }
	  f[i][0] = f[i][0] + delx*force;
	  f[i][1] = f[i][1] + dely*force;
	}
      }
    }
  }
	
}


/* Move the elk using Newtonian mechanics. */

void move_elk(n,x,v,f,radius)
  int    n;
  double x[MAX_ELK][2],v[MAX_ELK][2],f[MAX_ELK][2];
  double radius;
{

  int    i;
  double radsq,distsq,dt;

/* Simple 1st-order Euler integrator. */

  dt = 0.5;

  for (i = 0; i < n; i++) {
    x[i][0] = x[i][0] + dt*v[i][0] + 0.5*dt*dt*f[i][0];
    x[i][1] = x[i][1] + dt*v[i][1] + 0.5*dt*dt*f[i][1];
    v[i][0] = v[i][0] + dt*f[i][0];
    v[i][1] = v[i][1] + dt*f[i][1];
  }

/* Set velocity to 0.0 for elk that wander onto rancher's land. */

  radsq = radius * radius;

  for (i = 0; i < n; i++) {
    distsq = x[i][0]*x[i][0] + x[i][1]*x[i][1];
    if (distsq > radsq) {
      v[i][0] = 0.0;
      v[i][1] = 0.0;
    }
  }

}

/* Dump snapshot of coords to file. */

void dump(n,x,step,radius)
  int    n,step;
  double x[MAX_ELK][2];
  double radius;
{

  int    i;

  if (step == -1) {
    fclose (fp);
    return;
  }

  if (step == 0) {
    fp = fopen("dump.out","w");
    fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
    fprintf(fp,"%d\n",n);
    fprintf(fp,"ITEM: BOX BOUNDS\n");
    fprintf(fp,"%g %g\n",-2.0*radius,2.0*radius);
    fprintf(fp,"%g %g\n",-2.0*radius,2.0*radius);
  }
  
  fprintf(fp,"ITEM: TIMESTEP\n");
  fprintf(fp,"%d\n",step);
  fprintf(fp,"ITEM: ATOMS\n");
  for (i = 0; i < n; i++)
    fprintf(fp,"%d %d %g %g\n",i+1,1,x[i][0],x[i][1]);

}

/* Park/Miller RNG */
 
double rand1(iseed)
  int *iseed;
{

  double aa=16807.0;
  double mm=2147483647.0;
  double sseed;
  int    jseed;

  jseed = *iseed;
  sseed = jseed;
  jseed = aa*sseed/mm;
  sseed = aa*sseed - mm*jseed;
  jseed = sseed;
  
  *iseed = jseed;

  return(sseed/mm);

}