From 901b9f16494f37890be17ef4bb66e6efc6873340 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 10 Feb 2017 12:18:11 -0500
Subject: changed code to rely on jst

---
 src/spheres_poly/box.cpp | 1157 ----------------------------------------------
 1 file changed, 1157 deletions(-)
 delete mode 100644 src/spheres_poly/box.cpp

(limited to 'src/spheres_poly/box.cpp')

diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp
deleted file mode 100644
index 716090d..0000000
--- a/src/spheres_poly/box.cpp
+++ /dev/null
@@ -1,1157 +0,0 @@
-/* 
-   Updated July 24, 2009 to include hardwall boundary condition option, as 
-   well as polydisperse spheres.  
-
-   Packing of hard spheres via molecular dynamics
-   Developed by Monica Skoge, 2006, Princeton University
-   Contact: Monica Skoge (mskoge.princeton.edu) with questions
-   This code may be used, modified and distributed freely.
-   Please cite:
-   
-   "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
-   	M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
-	
-   if you use these codes.	
-*/
-
-#include "box.h"
-#include <iostream>
-#include <math.h>
-#include <stdlib.h>
-#include <time.h>
-#include <iomanip>
-
-//==============================================================
-//==============================================================
-//  Class Box: Fills box with hardspheres to given packing fraction
-//  and evolves spheres using molecular dynamics!
-//==============================================================
-//==============================================================
-
-
-//==============================================================
-// Constructor
-//==============================================================
-box::box(int N_i, double r_i, double growthrate_i, double maxpf_i, 
-	 double bidispersityratio_i, double bidispersityfraction_i, 
-	 double massratio_i, int hardwallBC_i):
-  r(r_i),          
-  N(N_i),
-  growthrate(growthrate_i),
-  h(N_i+1),
-  maxpf(maxpf_i),
-  bidispersityratio(bidispersityratio_i),
-  bidispersityfraction(bidispersityfraction_i),
-  massratio(massratio_i),
-  hardwallBC(hardwallBC_i)
-{
-  ngrids = Optimalngrids2(r);
-  cells.set_size(ngrids);
- 
-  s = new sphere[N];
-  binlist = new int[N];
-  x = new vector<DIM>[N];        
-  h.s = s;
-
-  gtime = 0.;
-  rtime = 0.;
-  ncollisions = 0;
-  ntransfers = 0;
-  nchecks = 0;
-  neventstot = 0;
-  ncycles = 0;
-  xmomentum = 0.; 
-  pressure = 0.;
-  collisionrate = 0.;
-
-  cells.set_size(ngrids);
-  cells.initialize(-1);      // initialize cells to -1
-  srandom(::time(0));        // initialize the random number generator
-  for (int i=0; i<N; i++)    // initialize binlist to -1
-    binlist[i] = -1;
-  
-  time(&start);
-}
-
-
-//==============================================================
-// Destructor
-//==============================================================
-box::~box() 
-{
-  delete[] s;
-  delete[] binlist;
-  delete[] x;
-}
-
-
-//==============================================================
-// ReadFile
-//==============================================================
-void box::ReadPositions(const char* filename)
-{  
-  // open file to read in arrays
-  std::ifstream infile(filename);
-  
-  infile.ignore(256, '\n');  // ignore the dim line
-  infile.ignore(256, '\n');  // ignore the #sphere 1 line
-  infile.ignore(256, '\n');  // ignore the #sphere line
-  infile.ignore(256, '\n');  // ignore the diameter line
-  infile.ignore(1000, '\n'); // ignore the 100 010 001 line
-  infile.ignore(256, '\n');  // ignore the T T T line
-
-  for (int i=0; i<N; i++)
-    {
-      infile >> s[i].r;      // read in radius    
-      infile >> s[i].gr;     // read in growth rate
-      infile >> s[i].m;      // read in mass
-      for (int k=0; k<DIM; k++)  
-	infile >> s[i].x[k]; // read in position 
-    }
-  infile.close();
-}
-
-
-//==============================================================
-// Recreates all N spheres at random positions
-//==============================================================
-void box::RecreateSpheres(const char* filename, double temp)
-{
-  ReadPositions(filename);  // reads in positions of spheres
-  VelocityGiver(temp);      // gives spheres initial velocities
-  AssignCells();            // assigns spheres to cells
-  SetInitialEvents();
-}
-
-
-//==============================================================
-// Creates all N spheres at random positions
-//==============================================================
-void box::CreateSpheres(double temp)
-{
-  int Ncurrent = 0;
-  for(int i=0; i<N; i++)
-    {
-      CreateSphere(Ncurrent);
-      Ncurrent++;
-    }
-  if (Ncurrent != N)
-    std::cout << "problem! only made " << Ncurrent << " out of " << N << " desired spheres" << std::endl;
-  
-  VelocityGiver(temp);
-  SetInitialEvents();
-}
-
-   
-//==============================================================
-// Creates a sphere of random radius r at a random unoccupied position 
-//==============================================================
-void box::CreateSphere(int Ncurrent)
-{
-  int keeper;    // boolean variable: 1 means ok, 0 means sphere already there
-  int counter = 0;   // counts how many times sphere already exists
-  vector<DIM> xrand;  // random new position vector
-  double d = 0.;
-  double radius;
-  double growth_rate;
-  double mass;
-  int species;
-  
-  if (Ncurrent < bidispersityfraction*N)
-    {
-      radius = r;
-      growth_rate = growthrate;
-      mass = 1.;
-      species = 1;
-    }
-  else
-    {
-      radius = r*bidispersityratio;
-      growth_rate = growthrate*bidispersityratio;
-      mass = massratio;
-      species = 2;
-    }
-  
-  while (counter<1000)
-    {
-      keeper = 1;
-      
-      for(int k=0; k<DIM; k++) 
-	xrand[k] = ((double)random()/(double)RAND_MAX)*SIZE;
-      
-      for (int i=0; i<Ncurrent; i++)  // check if overlapping other spheres
-	{
-	  d=0.;
-	  if (hardwallBC)
-	    {
-	      for (int k=0; k<DIM; k++) 
-		{     
-		  d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
-		}
-	    }
-	  else
-	    {
-	      for (int k=0; k<DIM; k++)
-		{
-		  if ((xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]) > SIZE*SIZE/4.)
-		    {
-		      if (xrand[k] > SIZE/2.)  // look at right image
-			d += (xrand[k] - (s[i].x[k]+SIZE))*
-			  (xrand[k] - (s[i].x[k]+SIZE));
-		      else                     // look at left image
-			d += (xrand[k] - (s[i].x[k]-SIZE))*
-			  (xrand[k] - (s[i].x[k]-SIZE));
-		    }
-		  else
-		    d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
-		}
-	    }
-	  if (d <= (radius + s[i].r)*(radius + s[i].r)) // overlapping!
-	    {
-	      keeper = 0;
-	      counter++;
-	      break;
-	    }
-	}
-      
-      if (hardwallBC)
-	{
-	  for (int k=0; k<DIM; k++)    // check if overlapping wall
-	    {     
-	      if ((xrand[k] <= radius) || (SIZE - xrand[k] <= radius)) // touching wall
-		{
-		  keeper = 0;
-		  counter++; 
-		  break;
-		}
-	    }
-	}
-      
-      if (keeper == 1)
-	break;   
-      
-      if (counter >= 1000)
-	{
-	  std::cout << "counter >= 1000" << std::endl;
-	  exit(-1);
-	}
-    }
-
-  // now convert xrand into index vector for cells
-  vector<DIM,int> cell;
-  cell = vector<DIM>::integer(xrand*((double)(ngrids))/SIZE);
-  
-  s[Ncurrent] = sphere(Ncurrent, xrand, cell, gtime, radius, growth_rate, 
-		       mass, species);
-  
-  //first check to see if entry at cell
-  if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
-    cells.get(cell) = Ncurrent;
-
-  else  // if no, add i to right place in binlist
-    {  
-      int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
-      int pointer = iterater;
-      while (iterater != -1)
-	{
-	  pointer = iterater;
-	  iterater = binlist[iterater];
-	}
-      binlist[pointer] = Ncurrent;
-    }
-
-}
-
-
-//==============================================================
-// Assign cells to spheres read in from existing configuration
-//==============================================================
-void box::AssignCells()
-{
-  for (int i=0; i<N; i++)
-    {
-      // now convert x into index vector for cells
-      vector<DIM,int> cell;
-      cell = vector<DIM>::integer(s[i].x*((double)(ngrids))/SIZE);
-      s[i].cell = cell;
-      
-      //first check to see if entry at cell
-      if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
-	cells.get(cell) = i;
-      
-      else  // if no, add i to right place in binlist
-	{  
-	  int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
-	  int pointer = iterater;
-	  while (iterater != -1)
-	    {
-	      pointer = iterater;
-	      iterater = binlist[iterater];
-	    }
-	  binlist[pointer] = i;
-	}
-    }
-}
-
-	
-//==============================================================
-// Velocity Giver, assigns initial velocities from Max/Boltz dist.
-//==============================================================
-void box::VelocityGiver(double T)
-{
-  for (int i=0; i<N; i++)
-    {
-      for (int k=0; k<DIM; k++)
-	{
-	  if (T==0.)
-	    s[i].v[k] = 0.;
-	  else
-	    s[i].v[k] = Velocity(T);
-	}
-    }
-}
-
-
-//==============================================================
-// Velocity, gives a single velocity from Max/Boltz dist.
-//==============================================================
-double box::Velocity(double T)
-{
-  double rand;                       // random number between -0.5 and 0.5
-  double sigmasquared = T;    // Assumes M = mass of sphere = 1
-  double sigma = sqrt(sigmasquared); // variance of Gaussian
-  double stepsize = 1000.;           // stepsize for discretization of integral
-  double vel = 0.0;                  // velocity
-  double dv=sigma/stepsize;
-  double p=0.0;
-  
-  rand = (double)random() / (double)RAND_MAX - 0.5;
-  if(rand < 0) 
-    {
-      rand = -rand;
-      dv = -dv;
-    }
-  
-  while(fabs(p) < rand) // integrate until the integral equals rand
-    {
-      p += dv * 0.39894228 * exp(-vel*vel/(2.*sigmasquared))/sigma;
-      vel += dv;
-    }
-  return vel;
-}
-
-
-//==============================================================
-// Finds next events for all spheres..do this once at beginning
-//==============================================================
-void box::SetInitialEvents()
-{
-  for (int i=0; i<N; i++)  // set all events to checks
-    {
-      event e(gtime, i, INF); 
-      s[i].nextevent = e;
-      h.insert(i);
-    }
-}
-
-
-//==============================================================
-// Finds next event for sphere i 
-//==============================================================
-event box::FindNextEvent(int i)
-{
-  double outgrowtime;
-  outgrowtime = (.5/ngrids - (s[i].r+s[i].gr*gtime))/s[i].gr + gtime;
-
-  event t = FindNextTransfer(i);
-  event c = FindNextCollision(i);
-
-  if ((outgrowtime < c.time)&&(outgrowtime < t.time)&&(ngrids>1))
-    {
-      event o = event(outgrowtime,i,INF-1);
-      return o;
-    }
-
-  if ((c.time < t.time)&&(c.j == INF)) // next event is check at DBL infinity
-    return c;
-  else if (c.time < t.time) // next event is collision!
-    {
-      CollisionChecker(c); 
-      return c; 
-    }
-  else // next event is transfer!
-    return t;  
-} 
-
-
-//==============================================================
-// Checks events of predicted collision partner to keep collisions
-// symmetric
-//==============================================================
-void box::CollisionChecker(event c)
-{
-  int i = c.i;
-  int j = c.j;
-  event cj(c.time,j,i,c.v*(-1));
-
-  // j should have NO event before collision with i!
-  if (!(c.time  < s[j].nextevent.time))
-    std::cout << i << " " <<  j << " error collchecker, s[j].nextevent.time= " << s[j].nextevent.time << " " << s[j].nextevent.j << ", c.time= " << c.time << std::endl;
-  
-  int k = s[j].nextevent.j; 
-  if ((k < N) && (k!=i)) // j's next event was collision so give k a check
-    s[k].nextevent.j = INF;
-  
-  // give collision cj to j
-  s[j].nextevent = cj;
-  h.upheap(h.index[j]);
-}
-
-
-//==============================================================
-// Find next transfer for sphere i 
-//==============================================================
-event box::FindNextTransfer(int i)
-{
-  double ttime = dblINF;  
-  int wallindex = INF;   // -(k+1) left wall, (k+1) right wall
-
-  vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
-  vector<DIM> vi = s[i].v;
-
-  for (int k=0; k<DIM; k++)
-    {
-      double newtime;
-      if (vi[k]==0.) 
-	newtime= dblINF;
-      else if (vi[k]>0)  // will hit right wall, need to check if last wall
-	{
-	  if ((hardwallBC)&&(s[i].cell[k] == ngrids - 1))
-	    newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
-		       - (xi[k]+s[i].r+s[i].gr*gtime))/(vi[k]+s[i].gr);
-	  else
-	    newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
-		       - xi[k])/(vi[k]);
-	  
-	  if (newtime<ttime)
-	    {
-	      wallindex = k+1;
-	      ttime = newtime;
-	    }
-	}
-      else if (vi[k]<0)  // will hit left wall
-	{
-	  if ((hardwallBC)&&(s[i].cell[k] == 0))
-	    newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) 
-		       - (xi[k]-(s[i].r+s[i].gr*gtime)))/(vi[k]-s[i].gr);
-	  else
-	    newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) 
-		       - xi[k])/(vi[k]);
-	  
-	  if (newtime<ttime)
-	    {
-	      wallindex = -(k+1);
-	      ttime = newtime;
-	    }
-	}
-    }
-  
-
-  if (ttime < 0)
-    ttime = 0;
-  // make the event and return it
-  event e = event(ttime+gtime,i,wallindex+DIM+N+1);
-  return e;
-}
-
-
-//==============================================================
-// Check all nearest neighbor cells for collision partners
-//==============================================================
-void box::ForAllNeighbors(int i, vector<DIM,int> vl, vector<DIM,int> vr,
-			  neighbor& operation)
-{
-  vector<DIM,int> cell = s[i].cell;
-
-  // now iterate through nearest neighbors
-  vector<DIM, int> offset;          // nonnegative neighbor offset
-  vector<DIM, int> pboffset;        // nearest image offset
-
-   vector<DIM,int> grid;
-
-   int ii ;
-
-   grid=vl;
-   while(1)
-   {
-     //if (vr[0] > 1)
-     //std::cout << grid << "..." << cell+grid << "\n";
-     for(int k=0; k<DIM; k++)
-     {
-        offset[k]=grid[k]+ngrids;  // do this so no negatives 	
-        if (cell[k]+grid[k]<0) //out of bounds to left
-          pboffset[k] = -1;
-	else if (cell[k]+grid[k]>=ngrids) // out of bounds to right
-	  pboffset[k] = 1;
-        else
-          pboffset[k] = 0;
-     }     
-     int j = cells.get((cell+offset)%ngrids);
-     while(j!=-1)
-       {
-	 operation.Operation(j,pboffset);
-	 j = binlist[j];
-       }
-
-     // A. Donev:     
-     // This code makes this loop dimension-independent
-     // It is basically a flattened-out loop nest of depth DIM
-     for(ii=0;ii<DIM;ii++)
-     {
-       grid[ii] += 1;
-       if(grid[ii]<=vr[ii]) break;
-       grid[ii]=vl[ii];
-     }
-     if(ii>=DIM) break;
-   }  
-}
-
-
-//==============================================================
-// PredictCollision
-//==============================================================
-void box::PredictCollision(int i, int j, vector<DIM, int> pboffset, 
-			     double& ctime, int& cpartner, 
-			     vector<DIM, int>& cpartnerpboffset)
-{
-  double ctimej;
-  
-  if (i!=j)
-    {	 
-      ctimej = CalculateCollision(i,j,pboffset.Double())+gtime;
-      
-      if (ctimej < gtime)
-	std::cout << "error in find collision ctimej < 0" << std::endl;
-      
-      if ((ctimej < ctime)&&(ctimej < s[j].nextevent.time))
-	{
-	  ctime = ctimej;
-	  cpartner = j;
-	  cpartnerpboffset = pboffset;
-	}	
-    }
-}
-
-
-//==============================================================
-// Find next collision
-//==============================================================
-event box::FindNextCollision(int i)
-{
-  collision cc(i, this);
-  
-  vector<DIM, int> vl, vr;
-
-  for (int k=0; k<DIM; k++)  // check all nearest neighbors
-    {
-      vl[k] = -1;
-      vr[k] = 1;
-    }
-  
-  ForAllNeighbors(i,vl,vr,cc);
-
-  event e;
-  if (cc.cpartner == i)  // found no collisions in neighboring cells
-    {
-      if (cc.ctime != dblINF)
-	std::cout << "ctime != dblINF" << std::endl;
-      e = event(dblINF,i,INF);  // give check at double INF
-    }
-  else
-    e = event(cc.ctime,i,cc.cpartner,cc.cpartnerpboffset);
-
-  return e;
-}
-
-
-//==============================================================
-// Calculates collision time between i and image of j using quadratic formula
-//==============================================================
-double box::CalculateCollision(int i, int j, vector<DIM> pboffset)
-{
-  if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12))
-    {
-      return dblINF;
-    }
-  else
-    {      
-      // calculate updated position and velocity of i and j
-      vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
-      vector<DIM> vi = s[i].v;
-      vector<DIM> xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime);
-      vector<DIM> vj = s[j].v;
-      
-      double r_now_i = s[i].r + gtime*s[i].gr;
-      double r_now_j = s[j].r + gtime*s[j].gr;
-      
-      double A,B,C;
-      A = vector<DIM>::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr);
-      B = vector<DIM>::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr);
-      C = vector<DIM>::norm_squared(xi - xj) - (r_now_i+r_now_j)*(r_now_i+r_now_j);
-
-      if (C < -1E-12*(r_now_i+r_now_j))
-	{
-	  std::cout << "error, " << i << " and " << j << 
-	    " are overlapping at time "<< gtime << std::cout;
-	  std::cout << "A, B, C = "  << A << " " << " " << B << 
-	    " " << " " << C <<  std::endl;
-	  if (CheckSphereDiameters()>0)
-	    std::cout << "a sphere has grown greater than unit cell" << 
-	      std::endl;
-	  else
-	    std::cout << "unknown error" << std::cout;
-	  exit(-1);
-	}
-      
-      return QuadraticFormula(A, B, C);
-    }  
-}
-
-
-//==============================================================
-// Quadratic Formula ax^2 + bx + c = 0
-//==============================================================
- double box::QuadraticFormula(double a, double b, double c)
-{
-  double x = dblINF;
-  double xpos;
-  double xneg;
-  double det = b*b - a*c;
-
-  if (c <= 0.)
-    {
-      if(b < 0.) // spheres already overlapping and approaching
-	{
-	  //std::cout << "spheres overlapping and approaching" << std::endl;
-	  //std::cout << "# events= " << neventstot << std::endl;
-	  x = 0.;	
-	}
-    }
-  else if (det > -10.*DBL_EPSILON)
-    {
-      if (det < 0.)  // determinant can be very small for double roots
-	det = 0.;    
-      if (b < 0.)
-	x = c/(-b + sqrt(det));
-      else if ((a < 0.)&&(b > 0.))
-	x = -(b + sqrt(det))/a;
-      else
-	x = dblINF;
-    }
-  return x;
-}
-
-
-//==============================================================
-// Returns first event
-//==============================================================
-void box::ProcessEvent()
-{  
-  neventstot++;
-  // Extract first event from heap
-  int i = h.extractmax();   
-  event e = s[i].nextevent; // current event
-  event f;                  // replacement event
-
-  if ((e.j>=0)&&(e.j<N))  // collision!
-    {
-      ncollisions++;
-      //std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
-      Collision(e);
-      f = FindNextEvent(i);
-      s[i].nextevent = f;
-      h.downheap(1);
-      if (f.time < e.time)
-	{
-	  std::cout << "error, replacing event with < time" << std::endl;
-	  exit(-1);
-	}
-      
-      // make sure collision was symmetric and give j a check
-      if ((s[e.j].nextevent.j != i)||(s[e.j].nextevent.time != gtime))
-	{
-	  std::cout << "error collisions not symmetric" << std::endl;
-	  std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
-	  std::cout << "but " << e.j << " thinks it has " << s[e.j].nextevent.j<< " "  << s[e.j].nextevent.time << std::endl;
-	  exit(-1);
-	}
-      else  // give j a check
-	s[e.j].nextevent.j = INF;
-    }
-  else if (e.j==INF)      // check!  
-    {
-      nchecks++;
-      //std::cout << "check for " << e.i << " at time " << e.time << std::endl;
-      f = FindNextEvent(i);
-      s[i].nextevent = f;
-      h.downheap(1);
-    }
-  else if (e.j==INF-1)      // sphere outgrowing unit cell, decrease ngrids!
-    {
-      gtime = e.time;
-      Synchronize(false);
-      ngrids = ngrids - 1;
-//      std::cout << "need to reduce ngrids to " << ngrids << std::endl;   
-      ChangeNgrids(ngrids);
-      h.downheap(1);
-    }
-  else                    // transfer!
-    {
-      ntransfers++;
-      //std::cout << "transfer for " << e.i << " at time " << e.time << std::endl;
-      Transfer(e);
-      f = FindNextEvent(i);
-      s[i].nextevent = f;
-      h.downheap(1);
-      //r = FindNextEvent(i, e.j-N-DIM-1);
-      //if (f.time <= e.time)
-      if (f.time < e.time)
-	{
-	  std::cout << "error after transfer, replacing new event with < time" << " " << std::endl;
-	  std::cout << std::setprecision(16) << "e.time= " << e.time << ", f.time= " << f.time << ", f.i= " << f.i << ", f.j= " << f.j << "e.i= " << e.i << ", e.j= " << e.j << std::endl;
-	  std::cout << std::setprecision(16) << "difference= " << e.time - f.time << std::endl;
-	  exit(-1);
-	}
-    }
-}
-
-
-//==============================================================
-// Processes a collision
-//=============================================================
-void box::Collision(event e)
-{
-  double ctime = e.time;
-  int i = e.i;
-  int j = e.j;
-  vector<DIM,int> v = e.v;  // virtual image
-  gtime = ctime;
-
-  // Update positions and cells of i and j to ctime
-  s[i].x += s[i].v*(gtime-s[i].lutime);
-  s[j].x += s[j].v*(gtime-s[j].lutime);
-
-  // Check to see if a diameter apart
-  double r_sum = s[i].r + s[j].r + (s[i].gr+s[j].gr)*gtime;
-  double distance = vector<DIM>::norm_squared(s[i].x - s[j].x- v.Double()*SIZE) - r_sum*r_sum;
-  if (distance*distance > 10.*DBL_EPSILON)
-    std::cout << "overlap " << distance << std::endl;
-
-  s[i].lutime = gtime;
-  s[j].lutime = gtime;
-  
-  vector<DIM,double> vipar;          // parallel comp. vi
-  vector<DIM,double> vjpar;          // parallel comp. vj
-  vector<DIM,double> viperp;         // perpendicular comp. vi
-  vector<DIM,double> vjperp;         // perpendicular comp. vj
-
-  // make unit vector out of displacement vector
-  vector<DIM,double> dhat;
-  dhat = s[i].x - s[j].x - v.Double()*SIZE;  // using image of j!!
-  double dhatmagnitude = sqrt(dhat.norm_squared());
-  dhat /= dhatmagnitude;
-
-  vipar = dhat*vector<DIM>::dot(s[i].v, dhat);
-  vjpar = dhat*vector<DIM>::dot(s[j].v, dhat);
-  viperp = s[i].v - vipar;
-  vjperp = s[j].v - vjpar;
-
-  s[i].v = vjpar + dhat*(s[i].gr+s[j].gr)*2 + viperp;
-  s[j].v = vipar - dhat*(s[i].gr+s[j].gr)*2 + vjperp;
-
-  // momentum exchange
-  double xvelocity;   // exchanged velocity
-  xvelocity = vector<DIM>::dot(s[i].v - s[j].v, dhat) - (s[i].gr+s[j].gr);
-  xmomentum += xvelocity*dhatmagnitude*s[i].m*s[j].m*2/(s[i].m+s[j].m);
-}
-
-
-//==============================================================
-// Transfer, takes care of boundary events too
-//=============================================================
-void box::Transfer(event e)
-{
-  gtime = e.time;
-  int i = e.i;
-  int j = e.j;
-  int k=0;           // dimension perpendicular to wall it crosses
- 
-  // update position and lutime (velocity doesn't change)
-  s[i].x += s[i].v*(gtime-s[i].lutime);
-  s[i].lutime = gtime;
-
-  vector<DIM,int> celli;  // new cell for i
-  celli = s[i].cell;  // this is not redundant
-  
-  // update cell
-  if (j>N+DIM+1)  // right wall
-    {
-      k = j-N-DIM-2;
-      celli[k] = s[i].cell[k] + 1;
-
-      if (hardwallBC)
-	{
-	  // if in right-most cell, reflect back
-	  if (s[i].cell[k] == ngrids - 1)
-	    s[i].v[k] = -s[i].v[k] - 4.*s[i].gr;
-	  else
-	    {
-	      if (ngrids>1)
-		UpdateCell(i, celli); 
-	    }
-	}
-      else
-	{
-	  // if in right-most cell, translate x and cell
-	  if (s[i].cell[k] == ngrids - 1)
-	    {
-	      s[i].x[k] -= SIZE;
-	      celli[k] -= ngrids;
-	    }
-	  if (ngrids>1)
-	    UpdateCell(i, celli); 
-	}
-    }
-  else if (j<N+DIM+1)  // left wall
-    {
-      k = -j+N+DIM;
-      celli[k] = s[i].cell[k] - 1;
-
-      if (hardwallBC)
-	{
-	  // if in left-most cell, reflect back
-	  if (s[i].cell[k] == 0)
-	    s[i].v[k] = -s[i].v[k] + 4.*s[i].gr;
-	  else
-	    {
-	      if (ngrids>1)
-		UpdateCell(i, celli); 
-	    }
-	}
-      else
-	{
-	  // if in left-most cell, translate x and cell
-	  if (s[i].cell[k] == 0)
-	    {
-	      s[i].x[k] += SIZE;
-	      celli[k] += ngrids;
-	    }
-	  if (ngrids>1)
-	    UpdateCell(i, celli); 
-	}
-    }
-  else
-    std::cout << "error in Transfer" << std::endl;    
-  
-}
-
-
-//==============================================================
-// Updates cell of a sphere to time
-//=============================================================
-void box::UpdateCell(int i, vector<DIM,int>& celli)
-{
-  if (celli == s[i].cell)
-    std::cout << "error in update cell..shouldn't be the same" << std::endl;
-  
-  // delete i from cell array at cell
-
-  if (cells.get(s[i].cell) == i) 
-    {
-      if (binlist[i] == -1)
-	cells.get(s[i].cell) = -1;
-      else
-	{
-	  cells.get(s[i].cell) = binlist[i];
-	  binlist[i] = -1;
-	}
-    }
-
-  else if (cells.get(s[i].cell) == -1)
-    {
-      std::cout << "error " << i << " not in claimed cell UpdateCell" << std::endl;
-      OutputCells();
-    }
-
-  else  // if no, find i in binlist
-    {  
-      int iterater = cells.get(s[i].cell);
-      int pointer = iterater;
-      while ((iterater != i)&&(iterater != -1))
-	{
-	  pointer = iterater;
-	  iterater = binlist[iterater];
-	}
-      if (iterater == -1)  // got to end of list without finding i
-	{
-	  std::cout << "problem " << i << " wasn't in claimed, cell iterater = -1" << std::endl;
-	  OutputCells();
-	}
-      else  // we found i!
-	{
-	  binlist[pointer] = binlist[i]; 
-	  binlist[i] = -1;
-	}	  
-    } 
-
-  // now add i to cell array at celli
-  s[i].cell = celli;
-  
-  //first check to see if entry at celli
-  if (cells.get(celli) == -1) //if yes, add i to cells gridfield
-    cells.get(celli) = i;
-  else  // if no, add i to right place in binlist
-    {
-      int iterater = cells.get(celli);  // now iterate through to end and add i
-      int pointer = iterater;
-      while (iterater != -1)  // find the end of the list
-	{
-	  pointer = iterater;
-	  iterater = binlist[iterater];
-	}
-      binlist[pointer] = i;
-      binlist[i] = -1; // redundant
-    }
-}
-
-
-//==============================================================
-// Output event heap...purely used for debugging
-//==============================================================
-void box::OutputEvents()
-{
-  h.print();
-}
-
-
-//==============================================================
-// Output positions of spheres and their cells...purely used for debugging
-//==============================================================
-void box::OutputCells()
-{
-  for (int i=0; i<N; i++)
-    std::cout << i << " " << s[i].x << " " << s[i].v << " " << s[i].cell << std::endl;
-}
-
-
-//==============================================================
-// Update positions...purely for graphical display
-//==============================================================
-void box::TrackPositions()
-{
-  for (int i=0; i<N; i++)
-    x[i] = s[i].x + s[i].v*(gtime-s[i].lutime);
-}
-
-
-//==============================================================
-// Computes the total energy
-//==============================================================
-double box::Energy()
-{
-  double E=0;
-  for (int i=0; i<N; i++)
-    E += 0.5*s[i].m*s[i].v.norm_squared();
-
-  return E/N;
-}
-
-
-//==============================================================
-// Calculates the packing fraction
-//==============================================================
-double box::PackingFraction()
-{
-  double rfactor = 0.;
-  for (int i=0; i<N; i++)
-    {
-      rfactor += pow(s[i].r + gtime*s[i].gr, DIM);
-    }
-  double v = (rfactor*pow(sqrt(PI), DIM))/(exp(lgamma(1.+((double)(DIM))/2.)));
-  return v/(pow(SIZE, DIM));
-}
-
-
-//==============================================================
-// Checks to make sure all sphere diameters are less than dimension
-// of unit cell
-//==============================================================
-int box::CheckSphereDiameters()
-{
-  int offender = 0;
-  for (int i=0; i<N; i++)
-    {
-      if (s[i].r*2 > 1./ngrids){
-	offender = i;
-	break;
-      }
-    }
-  return offender;
-}
-
-
-//==============================================================
-// Change ngrids
-//==============================================================
-void box::ChangeNgrids(int newngrids)
-{
-  cells.set_size(newngrids);
-  cells.initialize(-1);      // initialize cells to -1
-  for (int i=0; i<N; i++)    // initialize binlist to -1
-    binlist[i] = -1;
-  AssignCells();
-  for (int i=0; i<N; i++)
-    s[i].nextevent = event(0., i, INF); 
-  Process(N); 
-}	 
-
-
-//==============================================================
-// Calculates the optimal ngrids for the initial configuration
-// and assumes that ngrids gets updated (reduced) as the packing
-// proceeds
-//==============================================================
-int box::Optimalngrids2(double currentradius)
-{
-  return (int)(1./(2.*currentradius));
-}
-
-
-//==============================================================
-// Calculates the optimal ngrids, assuming ngrids is not updated
-// automatically and is very conservative
-//==============================================================
-int box::Optimalngrids()
-{
-  double maxr;
-
-  maxr = pow(exp(lgamma(1.+((double)(DIM))/2.))*maxpf/
-	     (N*(bidispersityfraction + (1.-bidispersityfraction)*
-		 pow(bidispersityratio, DIM))), 
-	     1./DIM)/sqrt(PI);
-
-  return (int)(1./(2.*maxr));
-}
-
-
-//==============================================================
-// Processes n events
-//==============================================================
-void box::Process(int n)
-{
-  double deltat = gtime;
-  for (int i=0; i<n; i++)
-    {
-      ProcessEvent();
-    }
-  pf = PackingFraction();   // packing fraction
-  deltat = gtime - deltat;
-  double oldenergy = energy;
-  energy = Energy();        // kinetic energy
-
-  energychange = ((oldenergy - energy)/oldenergy)*100; // percent change in energy
-
-  if (deltat != 0.)
-    {
-      pressure = 1+xmomentum/(2.*energy*N*deltat);
-      collisionrate = ((double)(ncollisions))/deltat;
-    }  
-
-  // reset to 0
-  ncollisions = 0;
-  ntransfers = 0;
-  nchecks = 0;
-  xmomentum = 0.;
-  ncycles++;
-}
-
-
-//==============================================================
-// Prints statistics for n events
-//==============================================================
-void box::PrintStatistics()
-{
-  std::cout << "packing fraction = " << pf << std::endl;
-  std::cout << "gtime = " << gtime << std::endl; 
-  std::cout << "total time = " << rtime+gtime << std::endl;
-  std::cout << "kinetic energy = " << energy << std::endl;
-  std::cout << "total # events = " << neventstot << std::endl;
-  std::cout << "# events = " << ncollisions+ntransfers+nchecks << ", # collisions = " << ncollisions << ", # transfers = " << ntransfers << ", # checks =" << nchecks << std::endl;
-  std::cout << "growthrate = " << growthrate << std::endl;
-  std::cout << "collisionrate = " << collisionrate << std::endl;
-  std::cout << "reduced pressure = " << pressure << std::endl;
-  std::cout << "-----------------" << std::endl;
-}
-
-
-//==============================================================
-// Updates spheres to gtime, synchronizes, and can change growth rate
-//==============================================================
-void box::Synchronize(bool rescale)
-{
-  double vavg = sqrt(2.*M*energy);
-
-  for (int i=0; i<N; i++)
-    {
-      s[i].x = s[i].x + s[i].v*(gtime-s[i].lutime);
-      s[i].nextevent.time -= gtime;
-
-      if (s[i].nextevent.time < 0.)
-	std::cout << "error, event times negative after synchronization" << std::endl;
-      if (rescale == true)   // give everyone checks
-	{
-	  s[i].nextevent = event(0., i, INF); 
-	  s[i].v /= vavg;
-	}
-	  
-      s[i].lutime = 0.;
-      s[i].r += gtime*s[i].gr;   
-    }
-
-  //r += gtime*growthrate;       // r defined at gtime = 0
-  rtime += gtime;
-  gtime = 0.;
-  
-  if (rescale == true)
-    Process(N);
-}
-
-
-//==============================================================
-// Run time
-//==============================================================
-void box::RunTime()
-{
-  time(&end);
-  std::cout << "run time = " << difftime(end, start) << std::endl;
-}
-
-
-//==============================================================
-// Write configuration
-//==============================================================
-void box::WriteConfiguration(double *data)
-{
-  int count1;
-
-  if (gtime != 0.)   // synchronize spheres if not currently synchronized
-    Synchronize(false);
-      
-  
- for (int i=0; i<N; i++)  // output radius and positions
-   {
-     for (int k=0; k<DIM; k++)
-			 data[DIM * i + k] = s[i].x[k];
-   }
- 
-}
-- 
cgit v1.2.3-70-g09d2