/********************************************************
% 
% 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.
% 
%-----

*/

#include <math.h>
#include <stdio.h>
#include "mpi.h"
#include "global_grid.h"
#include "sim_consts.h"
#include "updatePositions.h"

/* creates arrays for data that will be sent to each processor
*/ 
void updatePositions(int tsi, double dt){

  int i,cn,exptCntr[9]; 
  double fx_tot, fy_tot, dv_x, dv_y, dp_x, dp_y; 
  FILE *Fp; 

  /* open a file to save all position update calculations into ... 
   */
  sprintf(procInfo.pFilename,"ParSimResults/update_pos_proc_%d_ts_%d.dat",procInfo.my_rank,tsi);
  Fp=fopen(procInfo.pFilename,"w"); 

  /* exptCntr is used to keep track of which particle's exported force we are considering */ 
  for(i = 0; i < 9; i++) exptCntr[i]=0;

  /* 
     Calculate new velocities and positions for each body under "home" control
     Update the velocity with an increment dv_i = ( f_i / m_i ) dt
     """""" """ position """" "" """"""""" dp_i = ( v_i + dv_i / 2 ) dt
  */

  for(i = 0; i < procInfo.NbPP; i++){ /* for all of my "home" particles */ 

    /* the force internal to each box */ 
    fx_tot = procInfo.fx[i]; fy_tot = procInfo.fy[i];
    /*printf("fx_tot=%10.6f; fy_tot=%10.6f\n",fx_tot,fy_tot);*/

    /* accumulate any forces obtained if this particle was exported to any neighbor */
    for( cn=0; cn < 9; cn++ ){
      if(cn==4) continue;
      if( procInfo.pExportQ[i][cn] ){
	fx_tot += procInfo.forcesOnExportedPart[cn][2*exptCntr[cn]  ]; 
	fy_tot += procInfo.forcesOnExportedPart[cn][2*exptCntr[cn]+1]; 
	exptCntr[cn]=exptCntr[cn]+1; /* increment our particle counter */ 
      }
    }

    /* compute the amount to update the velocity by ... assuming unit masses for now */
    /*dv_x = ( fx_tot / procInfo.m[i] ) * dt; 
      dv_y = ( fy_tot / procInfo.m[i] ) * dt;*/
    dv_x = ( fx_tot ) * dt; 
    dv_y = ( fy_tot ) * dt; 

    /* save intermediates ... */ 
    fprintf(Fp,"%20.10f %20.10f %20.10f %20.10f %20.10f %20.10f ",
	    procInfo.px[i], procInfo.py[i], procInfo.vx[i],procInfo.vy[i],fx_tot,fy_tot);
    
    /* compute the amount to update the position by ... using midpoint rule */ 
    /*dp_x = ( procInfo.vx[i] + dv_x/2 ) * dt;
      dp_y = ( procInfo.vy[i] + dv_y/2 ) * dt;*/
    dp_x = ( procInfo.vx[i] + dv_x ) * dt;
    dp_y = ( procInfo.vy[i] + dv_y ) * dt;

    /* save intermediates ... */ 
    fprintf(Fp,"%20.10f %20.10f\n",
	    dp_x,dp_y);
    
    /* actually update my home particles */ 
    procInfo.vx[i] += dv_x; procInfo.vy[i] += dv_y; 
    procInfo.px[i] += dp_x; procInfo.py[i] += dp_y; 

    /* save intermediates ... */ 
    /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n",
      procInfo.px[i], procInfo.py[i], procInfo.vx[i],procInfo.vy[i],fx_tot,fy_tot,dt);*/
    
    procInfo.fx[i] = procInfo.fy[i] = 0.0; /* reset for force vector */ 
 
  } /* for each home particle */ 

  fclose(Fp); 
  
}