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

A serial implementation of a restricted range n-body calculation in
two dimensions.  This code is used to verify the correctness parallel
version of the midpoint method (as such it has some specific constants
hardwired in)

*/

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

/*#define _CONST_FORCE_ 1*/
#define _REAL_FORCE_ 1

/* The maximum number of bodies to simulate */ 
#define MAX_N_BODIES 1000000
/*#define gravG 6.67e-11*/
#define gravG 1.0
double DISTANCE_LOWER_BOUND=0.01;
/* the "direct" short range interaction radius (forces will be done with a direct O(N^2) method */ 
double DIRECT_SR_FORCE_RADIUS=12.0; 

/* 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 {
  double x,y;
} p[MAX_N_BODIES], v[MAX_N_BODIES], f[MAX_N_BODIES];

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

/* the timestep size */ 
double dt=1.0;

int NP, NbPP, Nts, tNP; 

/* function declarations */ 
void readParallelParticles(int,int); 
void calculateForces(int);
void moveBodies(int);
void saveBodies(int); 

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

  int gbi, ti, dum1, dum2, dum3, dum4, dum5, dum6, dum7, dum8, dum9; 
  double t, xtmp,ytmp,xdtmp,ydtmp;
  clock_t start, stop;

  /* Check input arguments
     Inputs are: 
       NP (number of processors this data was generated on)
       Nts (number of timesteps to run)             
  */
#ifdef __FALSE__
  if( argc < 3 | argc > 3 ){
    fprintf(stderr,"Incorrect number of input arguments ... exiting\n"); 
    exit(1);
  }
#endif
  /* Read in the number of processors this data was generated on */ 
  NP = atoi(argv[1]); 
  if( NP <= 0 ){
    fprintf(stderr,"Input argument *NP* must be positive\n");
    exit(1); 
  }
  /*printf("Number of Processors (NP)=%d\n",NP); */

  /* 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); 
  }
  /*printf("Number of timesteps (Nts)=%d\n",Nts); */

  start = clock();
  while (start == clock()) {
    ;
  }
  start = clock();

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

  stop = clock();
  /* output time */ 
  /*printf("Program execution time is given by = %10.6f\n", (double)(stop - start) / (double)CLOCKS_PER_SEC );*/
  {
    FILE *fp; 
    char output[100];
    sprintf(output,"SerSimResults/serial_timing_%d_%d.dat",tNP,Nts); 
    fp = fopen(output, "w");
    if( !fp ){ exit(-1); }
    fprintf(fp,"elapsed time = %10.6f seconds\n", (double)(stop - start) / (double)CLOCKS_PER_SEC );
    fclose(fp); 
  }

  return 0;
}

/* -- Read in from files the particle positions & velocities 
   -- Initialize forces & masses
*/ 
void readParallelParticles(int ti, int NP){

  int gbi,fi,bi;
  double xtmp,ytmp,xdtmp,ydtmp;
  char filename[80]; 
  FILE *Fp,*Fpb;
  
  gbi=0; tNP=0; 
  for( fi=0; fi < NP; fi++ ){

    /* 
       Read ASCII text files... 
     */
#ifdef __FALSE__
    /*sprintf(filename,"ParSimResults/sim_proc_rnk_%d.dat",fi);*/
    sprintf(filename,"ParSimResults/sim_proc_rnk_%d_ts_%d.dat",fi,ti);
    Fp=fopen(filename,"r"); 
    if( Fp == NULL ){
      fprintf(stderr,"Error opening file!\n"); exit(1);
    }
    
    /* read in the number of particle this file contains ... */ 
    fscanf(Fp,"%d\n",&NbPP); 
    printf("Number of bodies per processor (NbPP)=%d\n",NbPP);

    tNP += NbPP;  

    for( bi=0; bi < NbPP; bi++ ){
      fscanf(Fp,"%f %f %f %f %d %d %d %d %d %d %d %d %d", 
	&xtmp, &ytmp, &xdtmp, &ydtmp, &dum1, &dum2, &dum3, &dum4, &dum5, &dum6, &dum7, &dum8, &dum9 );
      /*fscanf(Fp,"%f %f %f %f\n", &xtmp, &ytmp, &xdtmp, &ydtmp);*/
      p[gbi].x=xtmp; p[gbi].y=ytmp; v[gbi].x=xdtmp; v[gbi].y=ydtmp;

      /*printf("data read = %10.6f %10.6f %10.6f %10.6f\n", (p[gbi].x), (p[gbi].y), (v[gbi].x), (v[gbi].y) );*/

      f[gbi].x = 0.0; f[gbi].y = 0.0; 
      m[gbi]=1.0; 

      gbi++; 
    } /* for bi loop */ 
    close(Fp); 
#endif
    /* 
       end read ASCII files... 
    */ 

    /* 
       Read Binary files... 
    */ 
    sprintf(filename,"ParSimResults/Binaries/sim_proc_rnk_%d_ts_%d.bin.dat",fi,ti);
    Fpb=fopen(filename,"rb"); 
    if( Fpb == NULL ){
      fprintf(stderr,"Error opening file!\n"); exit(1);
    }
    fread((char*)&NbPP, 1, sizeof(int), Fpb);
    /*printf("Number of bodies per processor (NbPP)=%d\n",NbPP);*/

    tNP += NbPP;

    for( bi=0; bi < NbPP; bi++ ){
      fread((char*)&xtmp,  1, sizeof(double), Fpb); 
      fread((char*)&ytmp,  1, sizeof(double), Fpb); 
      fread((char*)&xdtmp, 1, sizeof(double), Fpb); 
      fread((char*)&ydtmp, 1, sizeof(double), Fpb); 
      p[gbi].x=xtmp; p[gbi].y=ytmp; v[gbi].x=xdtmp; v[gbi].y=ydtmp;

      /*printf("data read = %10.6f %10.6f %10.6f %10.6f\n", (p[gbi].x), (p[gbi].y), (v[gbi].x), (v[gbi].y) );*/

      f[gbi].x = 0.0; f[gbi].y = 0.0; 
      m[gbi]=1.0; 

      gbi++; 
    } /* for bi loop */ 
    close(Fpb); 
    /* end read binary files... */ 

    /*printf("-----\n");*/
    
  }
}

/* Calculate the forces on every pair of bodies */ 
void calculateForces(int tsi){

  int i,j;
  double dx,dy,dist, f12, r12_x, r12_y; 
  char filename[80]; 
  FILE *Fp;

  sprintf(filename,"SerSimResults/force_calc_ts_%d.dat",tsi);
  Fp=fopen(filename,"w"); 

  for( i=0; i<tNP; i++ ){
    for( j=i+1; j<tNP; j++ ){
      
      dx = p[j].x - p[i].x; dy = p[j].y - p[i].y; 

      /*dist = sqrt( dx*dx + dy*dy );*/
      dist = sqrtf( dx*dx + dy*dy );
      /*printf("dist[%d,%d]=%6.3f\n",i,j,dist);*/

      if( dist >= DIRECT_SR_FORCE_RADIUS ){ /* no force calculate for particles spaced too far apart */ 
	continue; 
      }

      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);*/
      f12 = gravG / (dist*dist);
      
      /* The {\em unit} vector from body #1 to body #2 */
      r12_x = dx / dist;
      r12_y = dy / dist; 
      
#ifdef _CONST_FORCE_
      f[i].x += 1.0/3.0;
      f[i].y += 1.0/5.0;
      
      f[j].x += 1.0/3.0;
      f[j].y += 1.0/5.0;
#endif
#ifdef _REAL_FORCE_
      f[i].x +=  f12*r12_x; /* a direct update ... the force on body i due to body j */ 
      f[i].y +=  f12*r12_y; 
      
      f[j].x += -f12*r12_x; /* a symmetric update ... the force on body j due to body i */ 
      f[j].y += -f12*r12_y;
#endif

      /* save intermediates ... smallest x coordinate ordered first */ 
      if( p[i].x < p[j].x ){
	/*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
		p[i].x, p[i].y,
		p[j].x, p[j].y, 
		dx,dy,dist,f12,r12_x,r12_y,f12*r12_x,f12*r12_y);*/
	fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
		p[i].x, p[i].y,
		p[j].x, p[j].y, 
		f12*r12_x,f12*r12_y);
      }else{
	/*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
		p[j].x, p[j].y,
		p[i].x, p[i].y, 
		-dx,-dy,dist,f12,-r12_x,-r12_y,-f12*r12_x,-f12*r12_y);*/
	fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
		p[j].x, p[j].y,
		p[i].x, p[i].y, 
		-f12*r12_x,-f12*r12_y);
      }

    } /* for j */ 
  } /* for i */ 

  fclose(Fp); 

}

/* 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(int tsi){

  int i;
  double dv_x,dv_y,dp_x,dp_y;
  char filename[80]; 
  FILE *Fp; 

  /* open a file to save all position update calculations into ... 
   */
  sprintf(filename,"SerSimResults/move_bodies_ts_%d.dat",tsi);
  Fp=fopen(filename,"w"); 
  
  for( i=0; i<tNP; i++ ){
    /*dv_x = ( f[i].x / m[i] ) * dt; 
      dv_y = ( f[i].y / m[i] ) * dt; */
    dv_x = ( f[i].x ) * dt; 
    dv_y = ( f[i].y ) * dt; 

    /* save intermediates ... */ 
    fprintf(Fp,"%20.10f %20.10f %20.10f %20.10f %20.10f %20.10f ",
	    p[i].x, p[i].y, v[i].x, v[i].y, f[i].x, f[i].y);
    
    /*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;

    /* save intermediates ... */ 
    fprintf(Fp,"%20.10f %20.10f\n",
	    dp_x, dp_y);
    
    v[i].x += dv_x; v[i].y += dv_y; 
    p[i].x += dp_x; p[i].y += dp_y; 

    /* save intermediates ... */ 
    /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
      p[i].x, p[i].y, v[i].x, v[i].y, f[i].x, f[i].y, dt);*/

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

  fclose(Fp); 

}


/* 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 */ 
  char fn[100]; /* size of constant string + new line + max integer size in characters */ 
  FILE* fh; 
  
  sprintf(fn,"SerSimResults/seq_sim_ts_%d.dat",ti);
  fh = fopen( fn, "w" ); 
  for( i=0; i<tNP; i++ ){
    fprintf(fh, "%15.9f %15.9f %15.9f %15.9f\n", p[i].x,p[i].y,v[i].x,v[i].y);
  }
  fclose(fh); 
}