/********************************************************
% 
% Written by:
% -- 
% John L. Weatherwax                2006-05-12
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

A serial implementation of an n-body calculation in two dimensions.
The user is allowed to simulate the trajectories of Nb randomly placed
particles in the plane for Nts timesteps.

Note: a potential drawback of this code is that there is no adaptive
timestepping.  This was done so be able to compare the speedup for
parallel versions of this base algorithm.  Each would be playing on
somewhat of the same "field" ... each executing the same {\em number}
of timesteps.

*/

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

/* The maximum number of bodies to simulate */ 
#define MAX_N_BODIES 1000000
/*#define gravG 6.67e-11*/
#define gravG 1.0
/*#define dt 0.01*/
#define dt 0.001
#define DISTANCE_LOWER_BOUND 0.01

/* For this simple application we will code most of our variables globally */
/* p is the position of each body, 
   v is the velocity of each body, and 
   f is the force applied to each body
*/
struct point {
  float x,y;
} p[MAX_N_BODIES], v[MAX_N_BODIES], f[MAX_N_BODIES];

/* Each bodies mass */ 
float m[MAX_N_BODIES];

int Nb; 

/* function declarations */ 
void calculateForces();
void moveBodies();
void saveBodies(); 

extern int
main(int argc, char* argv[]){

  int Nts, bi, ti; 
  float t;

  /* Set evil seed (initial seed) */ 
  srand((unsigned)0);
  /* srand((unsigned)time(NULL)); */ 

  /* Check input arguments */
  if( argc < 3 | argc > 3 ){
    fprintf(stderr,"Incorrect number of input arguments ... exiting\n"); 
    exit(1);
  }
  /* Read in the requested number of particles */
  Nb = atoi(argv[1]);
  if( Nb <= 0 || Nb > MAX_N_BODIES ){
    fprintf(stderr,"Input argument *Nb* must be positive and less than %d\n",MAX_N_BODIES); 
    exit(1); 
  }
  /* Read in the length of time to run our simulation */ 
  /* tFinal = atof(argv[2]);
  if( tFinal <= 0 ){
    fprintf(stderr,"Input argument *tFinal* must be positive\n");
    exit(1); 
    }*/
  /* It seems that it would be easier to compare algorithms on the basis 
     of the number of timesteps take by each */ 
  Nts = atoi(argv[2]);
  if( Nts < 0 ){ 
    fprintf(stderr,"Input argument *Nts* must be positive\n"); 
    exit(1); 
  }

  /* -- Generate random particle positions & velocities in the square [-1,+1]x[-1,+1]
     -- Initialize forces & masses
  */ 
  for( bi=0; bi < Nb; bi++ ){
    p[bi].x = 2*((float) rand()/RAND_MAX)-1; p[bi].y = 2*((float) rand()/RAND_MAX)-1; 
    /*printf("p[%d].x = %f; p[%d].y = %f\n",bi,p[bi].x,bi,p[bi].y);*/ 
    v[bi].x = 2*((float) rand()/RAND_MAX)-1; v[bi].y = 2*((float) rand()/RAND_MAX)-1; 
    /*printf("v[%d].x = %f; v[%d].y = %f\n",bi,v[bi].x,bi,v[bi].y);*/ 
    f[bi].x = 0.0; f[bi].y = 0.0; 
    m[bi]=0.1; 
  }

  /* The main loop */ 
  t = 0.0;
  saveBodies(0);
  for( ti=1; ti <= Nts; ti++ ){
    printf("Integrating from %f to %f...\n",t,t+dt);
    calculateForces(); 
    moveBodies(); 
    saveBodies(ti); 
    t+=dt; 
  }

  return 0;
}

/* Calculate the forces on every pair of bodies */ 
void calculateForces(){
  int i,j;
  float dist, f12, r12_x, r12_y; 

  for( i=0; i<Nb; i++ ){
    for( j=i+1; j<Nb; j++ ){
      
      dist = sqrt( (p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) );
      /*printf("dist[%d,%d]=%6.3f\n",i,j,dist);*/

      if( dist < 0.01 ){
	/* printf("particles too close together...mollifying\n");  */
	dist+=DISTANCE_LOWER_BOUND;
      }

      /* The magnitude of the force */ 
      f12 = gravG * m[i] * m[j] / (dist*dist);
      
      /* The {\em unit} vector from body #1 to body #2 */
      r12_x = (p[j].x-p[i].x)/dist; 
      r12_y = (p[j].y-p[i].y)/dist;
      
      f[i].x +=  f12*r12_x; /* a direct update */ 
      f[i].y +=  f12*r12_y; 
      
      f[j].x += -f12*r12_x; /* a symmetric update */ 
      f[j].y += -f12*r12_y;
    }
  }
}

/* calculate new velocities and positions for each body
   Update the velocity with an increment dv_i = ( f_i / m_i ) dt
   """""" """ position """" "" """"""""" dp_i = ( v_i + dv_i / 2 ) dt
*/ 
void moveBodies(){

  float dv_x,dv_y;
  float dp_x,dp_y;
  int i; 
  
  for( i=0; i<Nb; i++ ){
    dv_x = ( f[i].x / m[i] ) * dt; 
    dv_y = ( f[i].y / m[i] ) * dt; 
    
    /*dp_x = ( v[i].x + dv_x/2 ) * dt;
      dp_y = ( v[i].y + dv_y/2 ) * dt;*/
    dp_x = ( v[i].x + dv_x ) * dt;
    dp_y = ( v[i].y + dv_y ) * dt;
    
    v[i].x += dv_x; v[i].y += dv_y; 
    p[i].x += dp_x; p[i].y += dp_y; 

    f[i].x = f[i].y = 0.0; /* reset for vector */ 
  }
}


/* Prints the point locations to a file for plotting purposes
*/ 
void saveBodies(int ti){

  int i; 
  char fn[18+1+10]; /* size of constant string + new line + max integer size in characters */ 
  FILE* fh; 
  
  sprintf(fn,"sequential_ts_%0004d.dat",ti);
  fh = fopen( fn, "w" ); 
  for( i=0; i<Nb; i++ ){
    fprintf(fh, "%10.6f     %10.6f     %10.6f     %10.6f\n", p[i].x,p[i].y,v[i].x,v[i].y);
  }
  fclose(fh); 
}