//Benjamin J. Hildebrand - Kenyon College '03 
//Written in the fall of 2001.  I continued the work of Dan Kiepfer's summerscience project.  
//This program is written to create a histogram of the differences in position of particles in Brownian Motion.  

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

//These two definitions must must agree with their definitions in bsumsci1.c.
#define NumParts 1 
#define NumTimeSteps 10

//Declaration of Functions
int main ();
void getdata ();
double deltas ();
void getdelts ();
int histogram ();
void gethist ();

//Declaration of global variables.
double data [NumParts * NumTimeSteps];
double delt [NumParts * (NumTimeSteps - 1)];
int hist [32];
FILE *fin, *fout, *fd, *fh, *fout2;
double delta;

//Main Program 
int main ()
{
  fout = fopen ("delts.dat", "w+");
  fout2 = fopen ("histdata.dat", "w+");
  getdata ();
  deltas ();
  fclose(fout);
  getdelts ();
  histogram ();
  fclose(fout2);
  gethist ();
}

//Obtaining input data from Dan Kiepfer's modified sumsci1.c (bsumsci1.c) output.

void getdata ()
{
  int r = 0; 
  fin = fopen ("data4.dat", "r"); //change file name here to input different position information of particles.

  while (fscanf (fin, "%lf", &data [r]) != EOF) 
    {
      printf ("%le\n", data [r]);
      r ++;
    }
}

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

double deltas ()
{
  int partnum = 0;
  int q = 0;
  int j,k,i;

  printf ("What timestep interval of delta would you like from the above data, (0,1,2,3,...),\n(where 0 is a difference of one time step)?  ");
  scanf ("%d", &j);
  printf ("The deltas have been written to the file, delts.dat.\n");
  
  for (k = NumParts + q; k <= NumParts * (NumTimeSteps - j) + 1; q ++)
    {        
      for (i = partnum; i <= NumParts * (NumTimeSteps - NumParts) + 1; i += NumParts)
	{
	  delta = data [k+j] - data [i];
	  printf ("%le\n", delta); //would display deltas in output.
	  fprintf (fout, "%le\n", delta); //writes deltas to delts.dat.
	  k += NumParts;
	}
      partnum ++;     
      if (partnum < NumParts)
	{
	  k = NumParts + q + 1;
	}
    }
}

//Reads in deltas from delts.dat and puts them into an array (delt).

void getdelts ()
{
  int r = 0; 
  fd = fopen ("delts.dat", "r");
 
  while (fscanf (fd, "%lf", &delt [r]) != EOF)
    {
      //printf ("%le\n", delt [r]);
      r ++;
    }
}

//Creates a histogram of the deltas.  

int histogram ()
{
  int u;
  int bin = 1;
  int counter = 0;
  double p;
  double lowerbound;
  double LowerBound;
  double upperbound;
  double UpperBound;
  double interval;
  double midpoint;
  double numstep;

  printf ("Input the Lower Boundary of the histogram:  ");
  scanf ("%le", &LowerBound);
  printf ("Input the Upper Boundary of the histogram:  ");
  scanf ("%le", &UpperBound);
  
  interval = (UpperBound - LowerBound) / 32;
  printf ("The Interval is: %le\n", interval);
  for (p = LowerBound; p <= UpperBound; p += interval)
    {
      upperbound = p + interval;
      printf ("Upper = %le ", upperbound);
      //fprintf (fout2, "%le ", upperbound);
      lowerbound = p;
      printf ("Lower = %le\n", lowerbound);
      //fprintf (fout2, "%le\n", lowerbound);
      midpoint = ((upperbound - lowerbound)/2) + lowerbound;
      fprintf (fout2, "%le ", midpoint);

      for (u = 0; u <= NumParts * (NumTimeSteps - 1); u ++)
        {
	  if (delt [u] < upperbound && delt [u] >= lowerbound)
	  {
	    //printf ("%d ", u);
	    counter ++;
	  }
	}
      printf ("Bin %d = %d\n", bin, counter);
      numstep = (double)counter / (NumParts * (NumTimeSteps - 1));
      fprintf (fout2, "%le\n", numstep);
      bin ++;
      counter = 0;
    }
}

void gethist ()
{
  int c = 0; 
  fh = fopen ("histdata.dat", "r");
 
  while (fscanf (fh, "%lf", &hist [c]) != EOF)
    {
      c ++;
    }
}