/*
reports force as a function of distance in fofrXXXX.dat.

The file sumsci4.c also reports the F(r) data, and you should use that file instead.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

#define ELECTRONCHARGE 1.6021892e-19 //charge on an electron in coulombs
#define EPSILON_NOUGHT 8.8541853e-12 //permittivity constant
#define BOLTZMANN 1.3806503e-23
#define KAPPA (1/161E-9) //equal to 1/screening_length
#define NUMPARTICLES 2 //must be 2 for this particlar simulation
#define NUMRUNS 1
#define MAXCHARS 20
#define OUTPUTACCURACY 12 //number of digits to output to files
//some default values
#define STARTDISTANCE 1.5e-6
#define STARTNUMSTEPS 1000
#define STARTTEMPERATURE 293.0
#define STARTVISCOSITY 0.0010019
#define STARTRADIUS 0.326e-6
#define STARTTIMESTEP 3.3e-2
#define STARTDIELECTRIC 80.0
#define STARTZ 1991
#define STARTTIMEBETWEENREPORTS 3.3e-2
 

// declare some structures
struct vector_t {
  double x, y, z; // the components in each of the three dimensions
};

struct particle_t {
  struct vector_t r; //the current position of the particle
  struct vector_t dr; //the deviation of the particle in a single timestep
  struct vector_t initial; //initial position of the particle
  double radius; //radius of the particle: we assume spherical particles
  double diffusivity; //the diffusivity of this particle: depends on some
  long Z; //charge on particle, in multiples of the elementary charge
  //of the world's variables
};

struct world_t {
  double viscosity; //viscosity of fluid in 10's of poise (SI units)
  double temperature; //temperature in kelvin
  double dielectric; //dielectric constant of fluid
};

//declare functions
int main();
int getLine(char input[MAXCHARS]);
void getParameters();
double drandUniform(double min, double max);
double drandGaussian(double stdDev);
void moveParticles(struct particle_t p[]);
void brownianMove(struct particle_t p[]);
void reportData(struct particle_t p[], long stepNo);
void reportParams(struct particle_t p[], struct world_t fluid);
void interactParticles(struct particle_t p[], struct world_t fluid);

//declare some external variables
FILE* fp; //output file
FILE* fp2; //temporary output file
struct particle_t p[NUMPARTICLES]; // the particles
struct world_t fluid; //some parameters of the fluid
long numSteps; //number of timesteps to take
double timeStep; //length of each timestep in seconds
char dataFile[20]; //the file to store the outputted data
char paramsFile[20]; //file to store the parameters of the system

//main program segment
int main() {

  //variable declarations
  extern struct particle_t p[NUMPARTICLES];
  extern struct world_t fluid;
  extern long numSteps;
  extern double timeStep;
  extern FILE* fp;
  extern FILE* fp2;
  extern char dataFile[20];
  extern char paramsFile[20];
  char fileSuffix[5];
  int i, j, k;
  time_t timetemp;
  long numStepsDivTen;
  int progressBar;
  char tempStr[30];
  
  //some initialization code
  timetemp = time(NULL);
  srand(timetemp);
  
  numSteps = 1000;
  timeStep = 3.3e-3;

  fluid.viscosity = 0.0010019; //viscosity in poise divided by 10: SI
  fluid.temperature = 293.0; //kelvin
  fluid.dielectric = 80.0; // about dielectric constant of water at 293K
  
  for (i=0; i<NUMPARTICLES; i++) {
    p[i].initial.x = p[i].initial.y = p[i].initial.z = 0.0;
    p[i].r = p[i].dr = p[i].initial;
    p[i].initial.x = p[i].r.x = i*1.5e-6;
    p[i].radius = 0.5e-6;
    p[i].diffusivity = BOLTZMANN * fluid.temperature/(6*M_PI*fluid.viscosity*p[i].radius);
    p[i].Z = 12000;
  }
  //if the user doesn't like the defaults given above, he can change them
  getParameters(); 

  //set up the progress reporting
  numStepsDivTen = numSteps / 10;
  numStepsDivTen = numStepsDivTen == 0 ? 1 : numStepsDivTen;
  progressBar = 0;
  printf("diffusivity of particles: %.4e\n", p[0].diffusivity);

  //set up the filenames to be used for output
  gcvt(timetemp % 10000, 4, fileSuffix);
  strcat(dataFile, "data");
  strcat(dataFile, fileSuffix);
  strcat(dataFile, ".dat");
  strcat(paramsFile, "params");
  strcat(paramsFile, fileSuffix);
  strcat(paramsFile, ".txt");
  fp = fopen(dataFile,"w");
  tempStr[0] = '\0';
  strcat(tempStr, "fofr");
  strcat(tempStr, fileSuffix);
  strcat(tempStr, ".dat");
  fp2 = fopen(tempStr,"w");

  //so the initial positions get reported
  reportData(p, 0);
  printf("number of particles: %d\n", NUMPARTICLES);
  printf("calculating:\n");

  //the main loop: each loop is one full simulation run
  for (k=0; k<NUMRUNS; k++) {
    
    //initializations that need to occur before each simulation
    for (i=0; i<NUMPARTICLES; i++) {
      p[i].r = p[i].initial;
    }


    //the timestep loop. each iteration corresponds to one timestep
    for (i=0; i<numSteps; i++) {
      //first apply some functions that will change the 'dr' 
      // of each particle
      //brownianMove(p);
      interactParticles(p, fluid);
      
      
      // once all functions have been applied to the particles,
      // update their positions
      moveParticles(p);
      
      //output information about the step
      reportData(p, i+1);
      
      //show a progress bar
      if ((i % numStepsDivTen) == 9) {
	printf("%d%%\n", progressBar);
	progressBar += 10;
      }
    } //end of a single simulation
  } //end of all simulations

  //shut down
  printf("100%%\n");
  fclose(fp);
  fclose(fp2);
  fp = fopen(paramsFile, "w");
  reportParams(p, fluid);
  printf("done.\nwritten data to file: %s\n", dataFile);
  printf("written parameter info to file: %s\n", paramsFile);
}

int getLine(char input[MAXCHARS]) {
  char c;
  int len = 0;
  while ((c = getchar()) != '\n' && len < MAXCHARS-1) {
    input[len] = c;
    len++;
  }
  input[len] = '\0';
  //printf("%sEND%d\n", input, len);
  return len;
}

void getParameters() {
  char tempStr[20];
  int i;
  extern struct particle_t p[NUMPARTICLES];
  extern struct world_t world;
  extern long numSteps;
  extern double timeStep;

  printf("You will now be prompted for some parameters of the system.\n");
  printf("If you don't wish to enter a certain parameter, then the\n");
  printf("default value in brackets will be used.\n\n");
  printf("Number of timesteps to take (%d): ", numSteps);
  if (getLine(tempStr))
    numSteps = atof(tempStr);
  printf("Length of each timestep in seconds (%.3e): ", timeStep);
  if (getLine(tempStr))
    timeStep = atof(tempStr);
  printf("Temperature of fluid in K (%.2fK): ", fluid.temperature);
  if (getLine(tempStr))
    fluid.temperature = atof(tempStr);
  printf("Viscosity of fluid in poise (%.5e poise): ", fluid.viscosity*10);
  if (getLine(tempStr))
    fluid.viscosity = atof(tempStr)/10; //convert to SI units
  printf("Dielectric constant of fluid (%.2f): ", fluid.dielectric);
  if (getLine(tempStr))
    fluid.dielectric = atof(tempStr);
  printf("Radius of colloids in meters (%.3em): ", p[0].radius);
  if (getLine(tempStr))
    for (i=0; i<NUMPARTICLES; i++) 
      p[i].radius = atof(tempStr);
  printf("Charge on each particle (%d electon charges): ", p[0].Z);
  if (getLine(tempStr))
    for (i=0; i<NUMPARTICLES; i++)
      p[i].Z = atof(tempStr);

  //in case any variables that affect diffusivity have been changed
  for (i=0; i<NUMPARTICLES; i++) 
    p[i].diffusivity = BOLTZMANN * fluid.temperature/(6*M_PI*fluid.viscosity*p[i].radius);
    

}

//this function moves all the particles to their positions
//for the next timestep.
void moveParticles(struct particle_t p[]){
  int i;
  for (i=0; i<NUMPARTICLES; i++) {
    p[i].r.x += p[i].dr.x;
    p[i].r.y += p[i].dr.y;
    p[i].r.z += p[i].dr.z;

    p[i].dr.x = p[i].dr.y = p[i].dr.z = 0.0;
  }
  return;
}

//this function applies some brownian randomness to each
//particle
void brownianMove(struct particle_t p[]) {
  extern double timeStep;  
  int i, j, k;
  double stdDev;
  double xtemp;
  
  //main brownian loop. brownianmove each particle
  for (i=0; i<NUMPARTICLES; i++) { 
    
    //randomly move this particle in each of the three dimensions
    stdDev = sqrt(2 * p[i].diffusivity * timeStep);
    p[i].dr.x += drandGaussian(stdDev);
    p[i].dr.y += drandGaussian(stdDev);
    p[i].dr.z += drandGaussian(stdDev);

  }
  return;
}

//returns a uniformly distributed random number between min and max
double drandUniform(double min, double max) {
  return rand()*(max-min)/RAND_MAX+min;
}

//returns a random number distributed gaussianly, with the given SD
double drandGaussian(double stdDev) {
  static int noExtra = 1;
  static double gset;
  double fac, r, v1, v2;

  if (noExtra) {
    do {
      v1 = drandUniform(-1.0, 1.0);
      v2 = drandUniform(-1.0, 1.0);
      r = v1*v1 + v2*v2;
    } while (r >= 1.0);
    fac = sqrt(-2.0*log(r)/r);
    gset = v1*fac;
    noExtra = 0;
    return v2*fac*stdDev;
  } else {
    noExtra = 1;
    return gset*stdDev;
  }
}


//this function outputs some data about each particle
//either to the screen or a file.
void reportData(struct particle_t p[], long stepNo) {
  extern FILE* fp;
  extern double timeStep;
  double distanceSq;
  double distance;
  char tempStr[OUTPUTACCURACY + 12];
  int i;
  double average;

  average = 0;
  gcvt(timeStep * stepNo, 9, tempStr);
  fputs(strcat(tempStr, "\t"), fp);

  //report the distance between two particles
  distance = sqrt((p[0].r.x-p[1].r.x)*(p[0].r.x-p[1].r.x) + (p[0].r.y-p[1].r.y)*(p[0].r.y-p[1].r.y) + (p[0].r.z-p[1].r.z)*(p[0].r.z-p[1].r.z));
  gcvt(distance, OUTPUTACCURACY, tempStr);
  fputs(strcat(tempStr, ""), fp);

  //report how far each particle has gone in total
  for (i=0; i<NUMPARTICLES; i++) {
    //distanceSq = p[i].r.x*p[i].r.x;// + p[i].r.y*p[i].r.y + p[i].r.z*p[i].r.z;
    //distance = sqrt(p[i].r.x*p[i].r.x + p[i].r.y*p[i].r.y + p[i].r.z*p[i].r.z);

    /* plot each particles x,y,z coords
    //once for each dimension    
    distance = p[i].r.x;
    gcvt(distance, OUTPUTACCURACY, tempStr);
    fputs(strcat(tempStr, "\t"), fp);
    distance = p[i].r.y;
    gcvt(distance, OUTPUTACCURACY, tempStr);
    fputs(strcat(tempStr, "\t"), fp);
    distance = p[i].r.z;
    gcvt(distance, OUTPUTACCURACY, tempStr);
    fputs(strcat(tempStr, "\t"), fp);
    */
  }
  fputs("\n", fp);
  //average = average / NUMPARTICLES;
  //gcvt(average, OUTPUTACCURACY, tempStr);
  //fputs(strcat(tempStr, "\n"), fp);  
}

//reports the parameters of the simulation to a file
void reportParams(struct particle_t p[], struct world_t fluid) {
  extern char paramsFile[];
  extern long numSteps;
  extern double timeStep;
  FILE* fp;
  int i;
  char outStr[80];

  fp = fopen(paramsFile, "w");
  
  fputs("Number of particles: ", fp);
  gcvt(NUMPARTICLES, 5, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Number of timesteps: ", fp);
  gcvt(numSteps, 5, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Length of each timestep (seconds): ", fp);
  gcvt(timeStep, 7, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Temperature of fluid (Kelvin): ", fp);
  gcvt(fluid.temperature, 6, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Viscosity of fluid (poise):", fp);
  gcvt(fluid.viscosity*10, 6, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Dielectric constant of fluid:", fp);
  gcvt(fluid.dielectric, 6, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fputs("Radius of particles (meters): ", fp);
  gcvt(p[0].radius, 6, outStr);
  fputs(strcat(outStr, "\n"), fp);
  
  fputs("Diffusivity of particles (meters^2/sec): ", fp);
  gcvt(p[0].diffusivity, 6, outStr);
  fputs(strcat(outStr, "\n"), fp);

  fclose(fp);
  
}

void interactParticles(struct particle_t p[], struct world_t fluid) {
  #define ELECTRONCHARGE 1.6021892e-19 //charge on an electron in coulombs
  #define EPSILON_NOUGHT 8.8541853e-12 //permittivity constant
  #define KAPPA 1/161e-9 //equal to 1/screening_length
  
  extern FILE* fp2;
  extern double timeStep;
  static double constantTerm[NUMPARTICLES]; //stores complicated calculations' results
  static int constantsCalculated = 0;
  double dr; //the distance the particle will be moved due to the force from another particle
  double f; //the force from one particle acting on another
  double r; //distance between the two particles
  double expTerm; //stores an exp() term that reoccurs - saves calculation time
  static double sixPiViscosity; //stores a constant
  double xDiff, yDiff, zDiff; // difference between two particles' displacement in each dimension
  int i, j, k; //some loop vars
  char tempStr[OUTPUTACCURACY + 12];
  
  if (!constantsCalculated) {
    for (i=0; i<NUMPARTICLES; i++) 
      constantTerm[i] = p[i].Z * ELECTRONCHARGE * exp(KAPPA*p[i].radius)/(1+KAPPA*p[i].radius) / sqrt(fluid.dielectric) / sqrt(4*M_PI*EPSILON_NOUGHT);
    sixPiViscosity = 6 * M_PI * fluid.viscosity;
    constantsCalculated = 1;
  }
  
  //run through each particle, moving it due to the influence from every other particle
  for (i=0; i<NUMPARTICLES; i++) {
    for (j=0; j<NUMPARTICLES; j++) {
      if (i != j) { //ensure the particle doesn't act on itself
	//move particle i due to the force from j
	xDiff = p[i].r.x-p[j].r.x;
	yDiff = p[i].r.y-p[j].r.y;
	zDiff = p[i].r.z-p[j].r.z;
	//printf("XPOS: %.20e\n",p[i].r.x);
	//printf("xdiff: %.20e\n",xDiff);
	r = sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff);
	//printf("total R: %.20e\n",r);
	expTerm = exp(-KAPPA*r);
	//printf("expterm: %.20e\n", expTerm);
	f = constantTerm[i] * constantTerm[j] * (KAPPA*expTerm/r + expTerm/(r*r)); 
	//printf("force: %.20e\n",f);
	dr = timeStep*f/(sixPiViscosity * p[i].radius);
	//printf("dr: %.20e\n",dr);
	//now add the displacements in each direction to the particle's dr
	p[i].dr.x += dr*xDiff/r;
	p[i].dr.y += dr*yDiff/r;
	p[i].dr.z += dr*zDiff/r;
	//printf("dx: %.20e\n",p[i].dr.x);
	//printf("dy: %.20e\n",p[i].dr.y);
	//printf("dz: %.20e\n",p[i].dr.z);
	gcvt(r, OUTPUTACCURACY, tempStr);
	fputs(strcat(tempStr, "\t"), fp2);
	gcvt(f, OUTPUTACCURACY, tempStr);
	fputs(strcat(tempStr, "\t"), fp2);
	//printf("%e, one: %e\n",constantTerm[i] * constantTerm[j], constantTerm[i]);
  
      } //end if
    } //end j loop
  } //end i loop
  fputs("\n", fp2);
}