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

*/

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

#define NUMPARTICLES 1
#define NUMSTEPS 80000 //number of timesteps to take
#define TIMESTEP 0.03333 //length of each timestep in seconds
#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, FILE *fpdataFile);
void reportParams (struct particle_t p[], struct world_t fluid, FILE *fpparams);
double deltas (FILE *fpdeltFile);

//declare some external variables
struct particle_t p [NUMPARTICLES]; // the particles
struct world_t fluid; //some parameters of the fluid
double data [(NUMPARTICLES + 1) * NUMSTEPS];
int filenum;
int filemin;
int filemax;
int l = 0;

//main program segment
int main (int argc, char *argv [])
{
 
  //variable declarations
  extern struct particle_t p [NUMPARTICLES];
  extern struct world_t fluid;
  FILE *fpdataFile, *fpdeltFile, *fpparams;
  double delta;
  char dataFile [20]; //the file to store the output data
  char paramsFile [20]; //file to store the parameters of the system
  char deltFile [20]; //file to store the output deltas
  char fileSuffix [5];
  int i, j, k;
  int r, ii;
  time_t timetemp;
  
  //some initialization code
  timetemp = time (NULL);
  srand (timetemp);

  //Instructions for running bsumsci2.c correctly.
  if (argc != 3)
    {
      printf ("Usage: ./bsumsci2 startfilenumber endfilenumber\n");
      exit (1);
    }
  sscanf (argv[1],"%d",&filemin);
  sscanf (argv[2],"%d",&filemax);

 //For loop that runs main a defined amount of times.
  for (filenum = filemin; filenum <= filemax; filenum++)
    {
      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 = 0.312e-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(); 
      //printf("diffusivity of particles: %.4e\n", p[0].diffusivity);

      //set up the filenames to be used for output
      gcvt (filenum, 4, fileSuffix);
      strcpy (dataFile, "data");
      strcat (dataFile, fileSuffix);
      strcat (dataFile, ".dat");
      strcpy (paramsFile, "params");
      strcat (paramsFile, fileSuffix);
      strcat (paramsFile, ".txt");
      fpdataFile = fopen (dataFile,"w+");

      strcpy (deltFile, "delt");
      strcat (deltFile, fileSuffix);
      strcat (deltFile, ".dat");  
      fpdeltFile = fopen (deltFile,"w+");

      printf ("Number of Particles: %d\n", NUMPARTICLES);
      printf ("Number of Time Steps: %d\n", NUMSTEPS);
     
      //so the initial positions get reported
      reportData (p, 0, fpdataFile);
      ii = 0;
      for (k=0; k<NUMPARTICLES; k++)
	{
	  data [k] = p[k].r.x;
	  //printf ("%le\n", data [ii]);//prints initial conditions in data []
	  ii++;
	}
      
      //the main loop. each iteration corresponds to one timestep
      for (j=0; j<NUMSTEPS; j++)
	{
	  //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, j+1, fpdataFile);
	  
	  for (r=0; r<NUMPARTICLES; r++)
	    {
	      ii = (j+1) * NUMPARTICLES + r;
	      data [ii] = p[r].r.x;
	      //printf ("%le\n", data [ii]); //prints data [] data.
	      ii++;
	    }
	}
      printf ("100%% done.");
     
      //Deltas takes the differeces in positions and then outputs that data
      //to deltXXXX.dat.
      deltas (fpdeltFile);

      //shut down
      fclose (fpdataFile);
      fclose (fpdeltFile);
      fpparams = fopen (paramsFile, "w");
      reportParams(p, fluid, fpparams);
      printf ("Data written to file: %s\n", dataFile);
      printf ("Parameter info written to file: %s\n", paramsFile);
      printf ("Deltas have been written to file: %s\n\n", deltFile);
      fclose (fpparams);
      
    }//end of filenum for loop
}//end of main

//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[]) {
  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, FILE *fpdataFile) 
{
  double distanceSq;
  double distance;
  char tempStr[OUTPUTACCURACY + 12];
  int i;

  //report how far each particle has gone in total
  for (i=0; i<NUMPARTICLES; i++) 
    {
      distance = 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(distance, OUTPUTACCURACY, tempStr);
      fputs(strcat(tempStr, "\n"), fpdataFile);
    }
}

//reports the parameters of the simulation to a file
void reportParams (struct particle_t p[], struct world_t fluid, FILE *fpparams) 
{
  char paramsFile[20];
  int i;
  char outStr[80];
  
  fputs("Number of particles: ", fpparams);
  gcvt(NUMPARTICLES, 5, outStr);
  fputs(strcat(outStr, "\n"), fpparams);

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

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

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

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

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

}

//Determining the change in position of each particle for every time step of a chosen difference in time intervals.

double deltas (FILE *fpdeltFile)
{
  int partnum = 0;
  int q = 0;
  //int l;
  int m,n;
  double delta;

  //printf ("\nWhat timestep interval of delta would you like, (0,1,2,3,...),\n(where 0 is a difference of one time step)?  ");
  //scanf ("%d", &l);
  printf ("\n%d\n", l);
      for (m = NUMPARTICLES + q; m <= NUMPARTICLES * (NUMSTEPS - l) + 1; q ++)
	{        
	  for (n = partnum; n <= NUMPARTICLES * (NUMSTEPS - NUMPARTICLES) + 1; n += NUMPARTICLES)
	    {
	      delta = data [m+l] - data [n];
	      //printf ("%le\n", delta); //displays deltas in output.
	      fprintf (fpdeltFile, "%le\n", delta); //writes deltas to deltsXXXX.dat.
	      m += NUMPARTICLES;
	    }
	  partnum ++;     
	  if (partnum < NUMPARTICLES)
	    {
	      m = NUMPARTICLES + q + 1;
	    }
	}
      if (l<7)
	{
	  l++;
	}
      else
	{
	  l = 0;
	}
}