/*
Runs many times to report many r(t+tau) vs. r(t) data points. tau is the timeBetweenReports
length of time, not the timeStep time. This data is outputted to dataXXXX.dat, where
XXXX is the last four digist of the time in seconds since Jan 1, 1970.

All output is in the form of a floating point number in ascii form, delimited
by tabs across rows. \n characters imply a new row.

This file also outputs F(r) data to a file fofrXXXX.dat.

A third output file, paramsXXXX.txt, contains the parameters of the run, so that
the simulation can be repeated if needed.

This file implements a local adaptive timestep algorithm. It can be switched off.

The particles' relevant information is stored in the p[] array. Their current
displacements in p[].r. Every timestep, the particles interact with each other
and have Brownian forces acting on them. These effect onto the p[].dr variable.
Then, only at the end of every timestep, is p[].dr added to p[].r. The alternative,
just changing p[].r when the particle needs to be moved in the middle of the
timestep, won't result in the conservation of momentum in some cases.

*/
#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 MAXCHARS 20
#define OUTPUTACCURACY 12 //number of digits to output to files
#define START_DISTANCE 1.5e-6
#define START_NUMRUNS 1000
#define START_NUMSTEPS 2000
#define START_TEMPERATURE 293.0
#define START_VISCOSITY 0.0010019
#define START_RADIUS 0.326e-6
#define START_TIMESTEP 3.3e-4
#define START_DIELECTRIC 80.0
#define START_Z 1991
#define START_TIMEBETWEENREPORTS 3.3e-2
#define STEPSBETWEENBROWNIAN 1

// 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
};

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();
double vectorLength(struct vector_t r);
double projectLength(struct vector_t r1, struct vector_t r2);
int getLine(char input[MAXCHARS]);
void getParameters();
double drandUniform(double min, double max);
double drandGaussian(double stdDev);
void moveParticles(struct particle_t p[]);
//implements the brownian motion of the particles
void brownianMove(struct particle_t p[]);
//reports the parameters of the simulation just run
void reportParams(struct particle_t p[], struct world_t fluid);
//handles the interactions of the particles
void interactParticles(struct particle_t p[], struct world_t fluid);

//declare some external variables
FILE* fp; //output file
FILE* fp2; //a second output file
struct particle_t p[NUMPARTICLES]; // the particles
struct world_t fluid; //some parameters of the fluid
long numRuns; //the number of times to run the simulation
long numSteps; //number of timesteps to take per run
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
struct vector_t thisR, prevR; //variables to report. reported in "interactParticles"
double timeBetweenReports=START_TIMEBETWEENREPORTS;

//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* fp;
  extern char dataFile[20];
  extern char paramsFile[20];
  char fileSuffix[5];
  int i, j, k;
  time_t timetemp;
  long numStepsDivTen;
  int progressBar;
  extern struct vector_t thisR, prevR;
  extern double timeBetweenReports;
  double timeElapsed; //time elapsed between reports
  char tempStr[OUTPUTACCURACY + 10];
  struct vector_t nullVector;
  long numFileLines = 0; //number of lines written to the data file

  //some initialization code
  timetemp = time(NULL);
  srand(timetemp);
  
  nullVector.x = nullVector.y = nullVector.z = 0;
  numRuns = START_NUMRUNS;
  numSteps = START_NUMSTEPS;
  timeStep = START_TIMESTEP;

  fluid.viscosity = START_VISCOSITY; //viscosity in poise divided by 10: SI
  fluid.temperature = START_TEMPERATURE; //kelvin
  fluid.dielectric = START_DIELECTRIC; // about dielectric constant of water at 293K
  
  //set the initial positions, and put the particles there
  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*START_DISTANCE;
    p[i].radius = START_RADIUS;
    p[i].diffusivity = BOLTZMANN * fluid.temperature/(6*M_PI*fluid.viscosity*p[i].radius);
    p[i].Z = START_Z;
  }

  //if the user doesn't like the defaults given above, he can change them
  getParameters(); 
  printf("diffusivity of particles: %.4e\n", p[0].diffusivity);

  //set up the progress reporting
  numStepsDivTen = numSteps / 10;
  numStepsDivTen = numStepsDivTen == 0 ? 1 : numStepsDivTen;
  progressBar = 0;

  //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");

  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;
    }

    // random initial distance
    //p[1].r.x = drandUniform(1.5e-6, 5e-6);
    
    //reset the distances between the two particles
    thisR = prevR = nullVector;
    
    //ensure the first positions are reported
    //reportings occur when timeElapsed > timeBetweenReports
    timeElapsed = timeStep/2;
    prevR = p[1].r; //the initial vector between the particles

    //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
      if (i % STEPSBETWEENBROWNIAN == STEPSBETWEENBROWNIAN-1) brownianMove(p);

      //this function also sets "thisR" to the vector between the particles
      interactParticles(p, fluid); 

      if (timeElapsed > timeBetweenReports) { //only report the data if the time has come
	gcvt(i*timeStep, OUTPUTACCURACY, tempStr);
	//fputs(strcat(tempStr, "\t"), fp);
	gcvt(vectorLength(prevR), OUTPUTACCURACY, tempStr);
	fputs(strcat(tempStr, "\t"), fp);
	gcvt(projectLength(thisR, prevR), OUTPUTACCURACY, tempStr);
	fputs(strcat(tempStr, "\n"), fp);
	timeElapsed = timeStep/2; //reset the time passed since a reporting to effectively 0
	prevR = thisR;
	numFileLines++;
      }
	
      
      // once all functions have been applied to the particles,
      // update their positions
      
      moveParticles(p);
      
      //count the time elapsed between reports
      timeElapsed += timeStep;

    } //end of a single simulation
    
    printf("%d of %d simulations done.\n", k, numRuns);
  } //end of all simulations

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

//return the length of the vector provided
double vectorLength(struct vector_t r) {
  return sqrt(r.x*r.x + r.y*r.y + r.z*r.z);
}

//return the length of the projection of r1 onto r2
double projectLength(struct vector_t r1, struct vector_t r2) {
  return (r1.x*r2.x + r1.y*r2.y + r1.z*r2.z)/vectorLength(r2);
}

//return the distance between two vectors
double vectorDistance(struct vector_t r1, struct vector_t r2) {
  double dx, dy, dz;
  dx = r1.x-r2.x;
  dy = r1.y-r2.y;
  dz = r1.z-r2.z;
  return sqrt(dx*dx + dy*dy + dz*dz);
}

//reads a line from input.
int getLine(char input[MAXCHARS]) {
  char c;
  int len = 0;
  while ((c = getchar()) != '\n' && len < MAXCHARS-1) {
    input[len] = c;
    len++;
  }
  input[len] = '\0';
  return len;
}

//get the simulation parameters from the user
void getParameters() {
  char tempStr[20];
  int i;
  extern struct particle_t p[NUMPARTICLES];
  extern struct world_t world;
  extern long numRuns;
  extern long numSteps;
  extern double timeStep;
  extern double timeBetwenReports;
  
  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 times to run the simulation (%d): ", numRuns);
  if (getLine(tempStr))
    numRuns = atof(tempStr);
  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("Time between reports in seconds (%.3e): ", timeBetweenReports);
  if (getLine(tempStr))
    timeBetweenReports = 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 * STEPSBETWEENBROWNIAN);
    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;
  }
}


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

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

  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("Length of reporting time (seconds): ", fp);
  gcvt(timeBetweenReports, 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);
  
}

//causes the particles to interact
void interactParticles(struct particle_t p[], struct world_t fluid) {

//maximum difference in % between force this timestep and force next timestep
//before adaptive timestepping occurs
#define FORCE_THRESHHOLD 0.05
  
  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
  double xProp, yProp, zProp; //proportion of movement in each dimension
  int i, j, k; //some loop vars
  char tempStr[OUTPUTACCURACY + 12];
  extern struct vector_t thisR, prevR;
  static double prevF=0;  
  double nextR; 
  double nextF;
  double adaptiveTimeStep;
  int adaptiveCount;
  double totalDr;
  int makeUpSteps;
  double projectedR; //r' projected onto r

  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;
	r = sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff);

	//the external distance - used for plotting r(t) vs r(t+t`)
	thisR.x = xDiff;
	thisR.y = yDiff;
	thisR.z = zDiff;

	xProp = xDiff/r;
	yProp = yDiff/r;
	zProp = zDiff/r;
	f = constantTerm[i] * constantTerm[j] * (KAPPA/r + 1/(r*r)) * exp(-KAPPA*r); 
	dr = timeStep*f/(sixPiViscosity * p[i].radius);
	nextR = r+dr;
	nextF = constantTerm[i] * constantTerm[j] * (KAPPA/nextR + 1/(nextR*nextR)) * exp(-KAPPA*nextR);
	adaptiveTimeStep = timeStep;
	adaptiveCount = 0;
	//first figure out the length of the adaptive timestep by repeatedly
	//halving it until the %difference in forces is less than FORCETHRESHHOLD
	while ((f - nextF)/f > FORCE_THRESHHOLD) {
	  adaptiveTimeStep /= 2;
	  adaptiveCount++;
	  dr = adaptiveTimeStep*f/(sixPiViscosity * p[i].radius);
	  nextR = r+dr;
	  nextF = constantTerm[i] * constantTerm[j] * (KAPPA/nextR + 1/(nextR*nextR)) * exp(-KAPPA*nextR);
	}
	//then, according to the length of the adaptiveTimeStep, run the force
	//simulation until the total time elapsed equals the original timeStep
	totalDr = 0;
	makeUpSteps = pow(2, adaptiveCount);

	//enable the following two lines to disable adaptive time-stepping
	adaptiveTimeStep = timeStep;
	makeUpSteps = 1;
	for (k=0; k < makeUpSteps; k++) {
	  f = constantTerm[i] * constantTerm[j] * (KAPPA/r + 1/(r*r)) * exp(-KAPPA*r);
	  dr = adaptiveTimeStep*f/(sixPiViscosity * p[i].radius);
	  r += 2*dr;
	  totalDr += dr;
	  //report f(r)
	  if (j==0) {
	    gcvt(r-2*totalDr, OUTPUTACCURACY, tempStr);
	    //fputs(strcat(tempStr, "\t"), fp2);
	    gcvt(f, OUTPUTACCURACY, tempStr);
	    //fputs(strcat(tempStr, "\n"), fp2);
	  }
 	}
	//now add the displacements in each direction to the particle's dr
	
	p[i].dr.x += totalDr*xProp;
	p[i].dr.y += totalDr*yProp;
	p[i].dr.z += totalDr*zProp;
	
	  
     } //end if
    } //end j loop

  } //end i loop


}