/*
This file outputs to a file, dataXXXX.dat, the average distance
squared particles are from their starting points over time. The data rows have
fields "time" and "avg r^2".

This is a very early version of the program, so if such data is required, sumsci4.c
should be modified to report such data instead.

*/

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

#define NUMPARTICLES 100
#define NUMRUNS 1
#define BOLTZMANN 1.3806503e-23
#define MAXCHARS 20
#define OUTPUTACCURACY 12 //number of digits to output to files


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

  double radius; //radius of the particle: we assume spherical particles
  double density; //density of the particle
  double mass; //mass of the particle
  double diffusivity; //the diffusivity of this particle: depends on some
  //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
};

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

//declare some external variables
FILE* fp; //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 char dataFile[20];
  extern char paramsFile[20];
  char fileSuffix[5];
  int i, j, k;
  time_t timetemp;
  long numStepsDivTen;
  int progressBar;

  //some initialization code
  timetemp = time(NULL);
  srand(timetemp);
  
  numSteps = 100;
  timeStep = 1e-6;

  fluid.viscosity = 0.0010019; //viscosity in poise divided by 10: SI
  fluid.temperature = 293.0; //kelvin
  
  for (i=0; i<NUMPARTICLES; i++) {
    p[i].r.x = p[i].r.y = p[i].r.z = 0;
    p[i].dr = p[i].r;
    p[i].radius = 1e-6;
    p[i].density = 0; //assumed to be same as water: not required
    p[i].mass = 0; //not required
    p[i].diffusivity = BOLTZMANN * fluid.temperature/(6*M_PI*fluid.viscosity*p[i].radius);
  }
  //if the user doesn't like the defaults given above, he can change them
  getParameters(); 
  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");

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

  //the main 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);



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

  //shut down
  printf("100%%\n");
  fclose(fp);
  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("Radius of colloids in meters (%.3em): ", p[0].radius);
  if (getLine(tempStr))
    for (i=0; i<NUMPARTICLES; i++) {
      p[i].radius = atof(tempStr);
      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;
  char tempStr[OUTPUTACCURACY + 12];
  int i;
  double average;

  average = 0;
  gcvt(timeStep * stepNo, 9, tempStr);
  fputs(strcat(tempStr, "\t"), 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;
    //gcvt(distanceSq, OUTPUTACCURACY, tempStr);
    //fputs(strcat(tempStr, "\t"), fp);
    average += distanceSq;
  }
  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("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);
  
}