#include <math.h>
#include <assert.h>

#include "LatticeEnvironment.h"

// The Lattice Environment is not instantiated initially
bool LatticeEnvironment::isInstantiated = false;

// Pointer to the Lattice Environment
LatticeEnvironment * LatticeEnvironment :: lEnv;

#define NUCLEUS_RADIUS (xDim / 20.0)

#define INACCESSIBLE_SMALL_PARTICLES
//#define RADIAL_MICROTUBULES
//#define RANDOM_MICROTUBULES
#define PLASMA_MEMBRANE     //Line only needed for signalling simulation

const static int CYTOSOL_ENVIRONMENT = 0;
static const int PLASMA_MEMBRANE_ENVIRONMENT = 1;
static const int NUCLEUS_ENVIRONMENT = 2;
static const int INACCESSIBLE_ENVIRONMENT = 3;

int NUCLEUS_VOLUME_TO_EXPOSE = 0;

float PI = 3.14159265;

//////////////////////////////////////////////////////////////////////
//  Constructor for the Lattice Environment
//
LatticeEnvironment::LatticeEnvironment(){

	stateMatrix = new Matrix(numLatticeSites, NUM_META);
	updateMatrix = new Matrix(numLatticeSites, NUM_META);
	statisticsMatrix = new Matrix(0,NUM_META);
	
	// Create the 3 Dimensional array of Lattice Sites
	buildLatticeSpace();

	// Build the list of neighbors for each lattice site
	buildNeighborList();

	countAndReportAccessibleNucleus();

	// Build Compatible list
	buildMoleculeCompatibleEnvTable();

	buildEnvironmentTable();
	
	populateStateMatrix();

	buildDiffusionRateTable();

	buildReactionRateTable();

	//buildDiffusionPatterns();       line not needed for signalling simulation
}



//////////////////////////////////////////////////////////////////////
//  Get instance of the Lattice Environment
//
LatticeEnvironment * LatticeEnvironment::getLatticeEnv(){

	// If Lattice Environment already instantiated, then return existing Environment
	// Else create new Lattice Environment and return the pointer
	if(!isInstantiated){

		// Create new Lattice Environment
		LatticeEnvironment * new_lEnv = new LatticeEnvironment();

		// Assign pointer
		lEnv = new_lEnv;
		
		// Lattice Environment exists
		isInstantiated = true;
	}

	// Return the pointer to the Lattice Environment
	return lEnv;
}



//////////////////////////////////////////////////////////////////////
//  Construct the Lattice Space 
//
void LatticeEnvironment :: buildLatticeSpace(){
	
	// make a 3D array of LatticeSites (Pointers)
	arr = new LatticeSite ***[xDim];

	// Step through each block of the array
	
	// Construct the X Dimension
	for (int x = 0 ; x < xDim ; x ++){

		arr[x] = new LatticeSite ** [yDim];
		
		// Construct the Y Dimension
		for (int y = 0 ; y < yDim ; y ++){

			arr[x][y] = new LatticeSite * [zDim];
			
			// Construct the Z Dimension
			for(int z = 0 ; z < zDim ; z ++){
			
				// Create a new LatticeSite at location (i,j,k)
				Point pt(x,y,z);
				
				// Create a Lattice Site at the specified Point
				arr[x][y][z] = new LatticeSite(pt);
			}
		}
	}
	// Construct the environment for the Lattice Space
	buildLatticeEnvironment();	
}



void LatticeEnvironment :: buildNeighborList(){

        // Step through the array and for each Lattice Site
        for(int x = 0 ; x < xDim ; x ++){
        for(int y = 0 ; y < yDim ; y ++){
        for(int z = 0 ; z < zDim ; z ++){

        LatticeSite * current = arr[x][y][z];

        for(int i = -1 ; i <= 1 ; i ++){
        for(int j = -1 ; j <= 1 ; j ++){
        for(int k = -1 ; k <= 1 ; k ++){
        	Point p((float)i,(float)j,(float)k);
                float mag = p.distFromOrigin();
                p.add(*(current->getLocation()));
                int degree = -1;
                        
		if(mag == (float)sqrt(1.0)) degree = 0;
                else if (mag == (float)sqrt(2.0)) degree = 1;
                else if (mag == (float)sqrt(3.0)) degree = 2;

                if(degree >= 0 ){
                if(withinCell(p)){
                	LatticeSite * neighbor = getLatticeSite(p);
	                current->neighbors[degree].push_back(neighbor);
        	}
		}

        }
        }
        }

	}
        }
        }

}

//////////////////////////////////////////////////////////////////////
//  Construct the environment for the Lattice space
//
void LatticeEnvironment :: buildLatticeEnvironment(){

	string line;
	string keyword;
	string delimiter = "\t";

	latticeMacroFeatures = NULL;

	struct LatticeMacroFeature * macroFeature;

	ifstream inFile (INPUT_FILE_ENV);

	while (!inFile.eof())
	{	
		getline(inFile,line,'\n');

		if(line.length() == 0 || (line[0] == '/' && line[1] == '/'))
			// a comment
			continue;

		int x1,y1,z1,x2,y2,z2,tmp, env;
		
		int keyPos = 0;

		macroFeature = new LatticeMacroFeature;

		// Get the type of construct
		keyPos = line.find(delimiter);
		keyword = line.substr(0,keyPos);

		if (strcmp( keyword.c_str(), "xplane" ) == 0) {
			// format: "xplane <x> <y1> <z1> <y2> <z2> <envType>"
			line = line.substr(keyPos+delimiter.length());

			// Get the fixed X coordinate of the plane
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			x2 = x1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get y1 coordinate
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get z1 coordinate
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y2 coordinate
			keyPos = line.find(delimiter);	
			y2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get z2 coordinate
			keyPos = line.find(delimiter);	
			z2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_XPLANE;
		}
		else if (strcmp( keyword.c_str(), "yplane" ) == 0) {
			// format: "yplane <y> <x1> <z1> <x2> <z2> <envType>"
			line = line.substr(keyPos+delimiter.length());

			// Get the fixed y coordinate of the plane
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			y2 = y1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get x1 coordinate
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get z1 coordinate
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get x2 coordinate
			keyPos = line.find(delimiter);	
			x2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get z2 coordinate
			keyPos = line.find(delimiter);	
			z2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_YPLANE;
		}
		else if (strcmp( keyword.c_str(), "zplane" ) == 0) {
			// format: "zplane <z> <x1> <y1> <x2> <y2> <envType>"
			line = line.substr(keyPos+delimiter.length());

			// Get the fixed Z coordinate of the plane
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			z2 = z1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get x1 coordinate
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y1 coordinate
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get x2 coordinate
			keyPos = line.find(delimiter);	
			x2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y2 coordinate
			keyPos = line.find(delimiter);	
			y2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_ZPLANE;
		}
		else if (strcmp( keyword.c_str(), "solid" ) == 0) {
			// format: "solid <x1> <y1> <z1> <x2> <y2> <z2> <envType>"
			line = line.substr(keyPos+delimiter.length());

			// Get x1 coordinate
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y1 coordinate
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get z1 coordinate
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get x2 coordinate
			keyPos = line.find(delimiter);	
			x2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y2 coordinate
			keyPos = line.find(delimiter);	
			y2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get z2 coordinate
			keyPos = line.find(delimiter);	
			z2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_SOLID;
		}
		else if (strcmp( keyword.c_str(), "hollow" ) == 0) {
			// format: "hollow <x1> <y1> <z1> <x2> <y2> <z2> <envType>"
			line = line.substr(keyPos+delimiter.length());

			// Get x1 coordinate
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y1 coordinate
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get z1 coordinate
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get x2 coordinate
			keyPos = line.find(delimiter);	
			x2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get y2 coordinate
			keyPos = line.find(delimiter);	
			y2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get z2 coordinate
			keyPos = line.find(delimiter);	
			z2 = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());

			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_HOLLOW;
		}
		else if (strcmp( keyword.c_str(), "point" ) == 0) {
			// Get x coordinate
			keyPos = line.find(delimiter);	
			x1 = atoi(line.substr(0,keyPos).c_str());
			x2 = x1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get Y position
			keyPos = line.find(delimiter);	
			y1 = atoi(line.substr(0,keyPos).c_str());
			y2 = y1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get Z position
			keyPos = line.find(delimiter);	
			z1 = atoi(line.substr(0,keyPos).c_str());
			z2 = z1;
			line = line.substr(keyPos+delimiter.length());
			
			// Get the environment
			env = atoi(line.substr(0).c_str());

			macroFeature->featureType = LATTICE_POINT;
		}
		else {
			cout << "Unrecognized environment specifier \"" << keyword.c_str();
			cout << "\".\n";
			exit(1);
		}

		if (x1 < 0) x1 = 0;
		if (x1 >= xDim) x1 = xDim - 1;
		if (x2 < 0) x2 = 0;
		if (x2 >= xDim) x2 = xDim - 1;
		if (x1 > x2) { tmp = x1; x1 = x2; x2 = tmp; }

		if (y1 < 0) y1 = 0;
		if (y1 >= yDim) y1 = yDim - 1;
		if (y2 < 0) y2 = 0;
		if (y2 >= yDim) y2 = yDim - 1;
		if (y1 > y2) { tmp = y1; y1 = y2; y2 = tmp; }

		if (z1 < 0) z1 = 0;
		if (z1 >= zDim) z1 = zDim - 1;
		if (z2 < 0) z2 = 0;
		if (z2 >= zDim) z2 = zDim - 1;
		if (z1 > z2) { tmp = z1; z1 = z2; z2 = tmp; }

		int x, y, z;

		if (macroFeature->featureType == LATTICE_HOLLOW) {
			for (x = x1; x <= x2; x++)
				for (z = z1; z <= z2; z++)
					arr[x][y1][z]->setEnvironment(env);
			for (x = x1; x <= x2; x++)
				for (z = z1; z <= z2; z++)
					arr[x][y2][z]->setEnvironment(env);
			for (x = x1; x <= x2; x++)
				for (y = y1; y <= y2; y++)
					arr[x][y][z1]->setEnvironment(env);
			for (x = x1; x <= x2; x++)
				for (y = y1; y <= y2; y++)
					arr[x][y][z2]->setEnvironment(env);
			for (y = y1; y <= y2; y++)
				for (z = z1; z <= z2; z++)
					arr[x1][y][z]->setEnvironment(env);
			for (y = y1; y <= y2; y++)
				for (z = z1; z <= z2; z++)
					arr[x2][y][z]->setEnvironment(env);
		}
		else {
			for (int x = x1; x <= x2; x++)
				for (int y = y1; y <= y2; y++)
					for (int z = z1; z <= z2; z++)
						arr[x][y][z]->setEnvironment(env);
		}

		macroFeature->x1      = x1;
		macroFeature->x2      = x2;
		macroFeature->y1      = y1;
		macroFeature->y2      = y2;
		macroFeature->z1      = z1;
		macroFeature->z2      = z2;
		macroFeature->envType = env;

		macroFeature->next    = latticeMacroFeatures;
		latticeMacroFeatures  = macroFeature;
	}
	
	inFile.close();

	buildNucleus();

	NUCLEUS_VOLUME_TO_EXPOSE = getI("NUCLEUS_VOLUME_TO_EXPOSE");

	bool *** accessible_nucleus_lattice;

	accessible_nucleus_lattice = (bool ***) malloc (xDim * sizeof(bool **));
	int ix, iy, iz;

	for (ix = 0; ix < xDim; ix ++)
	{
	  accessible_nucleus_lattice[ix] = (bool **) malloc (yDim * sizeof (bool *));
	  for (iy = 0; iy < yDim; iy ++)
	  {
	    accessible_nucleus_lattice[ix][iy] = (bool *) malloc (zDim * sizeof (bool));
	    for (iz = 0; iz < zDim; iz ++)
	    {
	      accessible_nucleus_lattice[ix][iy][iz] = false;
	    }
	  }
	}

	int accessible_nucleus;

	accessible_nucleus = countAndRecordAccessibleNucleus(NULL, accessible_nucleus_lattice);

#ifdef INACCESSIBLE_SMALL_PARTICLES
	buildInaccessibleSpace();
#else
	buildMicrotubules(accessible_nucleus_lattice, accessible_nucleus - NUCLEUS_VOLUME_TO_EXPOSE);
#endif /* INACCESSIBLE_SMALL_PARTICLES */

#ifdef PLASMA_MEMBRANE
	buildPlasmaMembrane();
#endif /* PLASMA_MEMBRANE */

	for (ix = 0; ix < xDim; ix ++)
	{
	  for (iy = 0; iy < yDim; iy ++)
	  {
	    free (accessible_nucleus_lattice[ix][iy]);
	  }
	  free (accessible_nucleus_lattice[ix]);
	}
	free (accessible_nucleus_lattice);
	accessible_nucleus_lattice = NULL;
}

void LatticeEnvironment :: buildNucleus(){
  double center_x, center_y, center_z;
  
  center_x = xDim / 2.0;
  center_y = yDim / 2.0;
  center_z = zDim / 2.0;

  for (int x = 0; x <= xDim; x++)
    for (int y = 0; y <= yDim; y++)
      for (int z = 0; z <= zDim; z++)
      {
	double distance_from_center;
	
	distance_from_center = sqrt((center_x - x) * (center_x - x) + (center_y - y) * (center_y - y) + (center_z - z) * (center_z - z));

	if (distance_from_center <= NUCLEUS_RADIUS)
	{
	  arr[x][y][z]->setEnvironment(NUCLEUS_ENVIRONMENT);
	}

      }
}

void LatticeEnvironment :: buildInaccessibleSpace(){
  static const int INACCESSIBLE_SPACE_PERCENT = getI("INACCESSIBLE_SPACE_PERCENT");

  int num_inaccessible_cells = INACCESSIBLE_SPACE_PERCENT * xDim * yDim * zDim / 100;

  int c = 0;

  while (c < num_inaccessible_cells)
  {
  
    int x, y, z;
    x = RandomNG::randInt(0, xDim - 1);
    y = RandomNG::randInt(0, yDim - 1);
    z = RandomNG::randInt(0, zDim - 1);

    if (arr[x][y][z]->getEnvironment() == 0)
    {
      arr[x][y][z]->setEnvironment(INACCESSIBLE_ENVIRONMENT);
      c++;
    }
  }
}

void LatticeEnvironment :: buildMicrotubules(bool *** accessible_nucleus_lattice, int max_nucleus_to_block){
  static const int INACCESSIBLE_SPACE_PERCENT = getI("INACCESSIBLE_SPACE_PERCENT");

  int num_inaccessible_cells = INACCESSIBLE_SPACE_PERCENT * xDim * yDim * zDim / 100;

  int cells_blocked = 0;

#ifdef RADIAL_MICROTUBULES
  while (cells_blocked < num_inaccessible_cells)
  {

    Point cell_membrane_end;
    cell_membrane_end = pickRandomPointOnCellMembrane();

    cells_blocked += buildRadialCylinder(cell_membrane_end,
					 NUCLEUS_RADIUS,
					 1, INACCESSIBLE_ENVIRONMENT, num_inaccessible_cells - cells_blocked, accessible_nucleus_lattice, &max_nucleus_to_block);
  }
#endif /* RADIAL_MICROTUBULES */

#ifdef RANDOM_MICROTUBULES
  while (cells_blocked < num_inaccessible_cells)
  {
    Point end1;
    end1.x = RandomNG::randInt(1, xDim - 1);
    end1.y = RandomNG::randInt(1, yDim - 1);
    end1.z = RandomNG::randInt(1, zDim - 1);

    Point end2;
    end2.x = RandomNG::randInt(1, xDim - 1);
    end2.y = RandomNG::randInt(1, yDim - 1);
    end2.z = RandomNG::randInt(1, zDim - 1);

    cells_blocked += buildCylinder(end1, end2,
				   1, INACCESSIBLE_ENVIRONMENT, num_inaccessible_cells - cells_blocked, accessible_nucleus_lattice, &max_nucleus_to_block);
  }
#endif /* RANDOM_MICROTUBULES */
}

int LatticeEnvironment :: buildRadialCylinder(Point cell_membrane_end,
			 double endpoint_distance_from_center,
			 int radius, int environment, int max_cells_to_block, bool *** accessible_nucleus_lattice, int * max_nucleus_to_block_ptr)
{
  Point nucleus_end;

#if 0
  double unit_x, unit_y, unit_z;
  double center_x, center_y, center_z;
  
  center_x = xDim / 2.0;
  center_y = yDim / 2.0;
  center_z = zDim / 2.0;

  double distance_from_center;

  distance_from_center = sqrt(
			      (center_x - cell_membrane_end.x) * (center_x - cell_membrane_end.x)
			      + (center_y - cell_membrane_end.y) * (center_y - cell_membrane_end.y)
			      + (center_z - cell_membrane_end.z) * (center_z - cell_membrane_end.z));

  unit_x = (center_x - cell_membrane_end.x) / distance_from_center;
  unit_y = (center_y - cell_membrane_end.y) / distance_from_center;
  unit_z = (center_z - cell_membrane_end.z) / distance_from_center;

  double length;

  length = distance_from_center - endpoint_distance_from_center;
  
  nucleus_end.x = (int) (cell_membrane_end.x + unit_x * length);
  nucleus_end.y = (int) (cell_membrane_end.y + unit_y * length);
  nucleus_end.z = (int) (cell_membrane_end.z + unit_z * length);
#endif

  nucleus_end = pickRandomPointOnNucleus(endpoint_distance_from_center);

  return buildCylinder(cell_membrane_end, nucleus_end,
		       radius, environment, max_cells_to_block, accessible_nucleus_lattice, max_nucleus_to_block_ptr);
}

int LatticeEnvironment :: buildCylinder(Point end1, Point end2,
					int radius, int environment, int max_cells_to_block, bool *** accessible_nucleus_lattice, int * max_nucleus_to_block_ptr)
{
  int cells_blocked;

  cells_blocked = buildCylinder_or_check(end1, end2, radius, environment, max_cells_to_block, false, accessible_nucleus_lattice, max_nucleus_to_block_ptr);

  if (cells_blocked < max_cells_to_block && cells_blocked != -1)
  {
    return buildCylinder_or_check(end1, end2, radius, environment, max_cells_to_block, true, accessible_nucleus_lattice, max_nucleus_to_block_ptr);
  }

  return cells_blocked;
}

int LatticeEnvironment :: buildCylinder_or_check(Point end1, Point end2,
						 int radius, int environment, int max_cells_to_block,
						 bool build, bool *** accessible_nucleus_lattice, int * max_nucleus_to_block_ptr)
{
  int * blocked_nucleus_cells_x;
  int * blocked_nucleus_cells_y;
  int * blocked_nucleus_cells_z;

  blocked_nucleus_cells_x = (int *) malloc (sizeof (int) * (xDim * yDim * zDim));
  blocked_nucleus_cells_y = (int *) malloc (sizeof (int) * (xDim * yDim * zDim));
  blocked_nucleus_cells_z = (int *) malloc (sizeof (int) * (xDim * yDim * zDim));

  int cells_blocked = 0;
  int nucleus_blocked = 0;

  // Here's an inefficient attempt at building a solid, round-capped cylinder.
  // Microtubules are actually hollow, uncapped cylinders.

  double length;

  length = sqrt ((double) ((end2.x - end1.x) * (end2.x - end1.x) + (end2.y - end1.y) * (end2.y - end1.y) + (end2.z - end1.z) * (end2.z - end1.z)));

  double unit_x, unit_y, unit_z;

  unit_x = (end2.x - end1.x) / length;
  unit_y = (end2.y - end1.y) / length;
  unit_z = (end2.z - end1.z) / length;

  int units_along_axis;

  for (units_along_axis = 0; units_along_axis < length; units_along_axis ++)
  {
    double center_x, center_y, center_z;

    center_x = end1.x + unit_x * units_along_axis;
    center_y = end1.y + unit_y * units_along_axis;
    center_z = end1.z + unit_z * units_along_axis;

    int x, y, z;

    for (x = (int) center_x - radius - 1; x < (int) center_x + radius + 1; x ++)
      for (y = (int) center_y - radius - 1; y < (int) center_y + radius + 1; y ++)
	for (z = (int) center_z - radius - 1; z < (int) center_z + radius + 1; z ++)
	{
	  double distance_to_center;

	  distance_to_center = sqrt ((double) ((x - center_x) * (x - center_x) + (y - center_y) * (y - center_y) + (z - center_z) * (z - center_z)));
	  
	  if (distance_to_center <= radius)
	  {
	    if (x >= 0 && x < xDim &&
		y >= 0 && y < yDim &&
		z >= 0 && z < zDim)
	    {
	      if (cells_blocked < max_cells_to_block)
	      {
		if (arr[x][y][z]->getEnvironment() == CYTOSOL_ENVIRONMENT)
		{
		  if (build)
		  {
		    arr[x][y][z]->setEnvironment(INACCESSIBLE_ENVIRONMENT);
		  }

		  cells_blocked ++;
		  if (accessible_nucleus_lattice[x][y][z])
		  {
		    blocked_nucleus_cells_x[nucleus_blocked] = x;
		    blocked_nucleus_cells_y[nucleus_blocked] = y;
		    blocked_nucleus_cells_z[nucleus_blocked] = z;
		    nucleus_blocked ++;
		    accessible_nucleus_lattice[x][y][z] = false;
		  }
		}
	      }
	    }
	  }
	}
  }

  if (!build)
  {
    int i;
    Point i_point;
    for (i = 0; i < nucleus_blocked; i++)
    {
      int x, y, z;
      x = blocked_nucleus_cells_x[i];
      y = blocked_nucleus_cells_y[i];
      z = blocked_nucleus_cells_z[i];

      accessible_nucleus_lattice[x][y][z] = true;
    }
  }

  free (blocked_nucleus_cells_x);
  free (blocked_nucleus_cells_y);
  free (blocked_nucleus_cells_z);

  if (nucleus_blocked > (*max_nucleus_to_block_ptr))
  {
    return -1;
  }
  
  if (build)
  {
    (*max_nucleus_to_block_ptr) -= nucleus_blocked;
  }

  return cells_blocked;
}

void LatticeEnvironment :: buildPlasmaMembrane(){
  int x, y, z;

  for (x = 0; x < xDim; x ++)
    for (y = 0; y < yDim; y ++)
      for (z = 0; z < zDim; z ++)
      {
	if (x == 0 || x == xDim - 1 || y == 0 || y == yDim - 1 || z == 0 || z == zDim - 1)
	{
	  arr[x][y][z]->setEnvironment(PLASMA_MEMBRANE_ENVIRONMENT);
	}
      }

}

Point LatticeEnvironment :: pickRandomPointOnNucleus(double nucleus_radius)
{
  double center_x, center_y, center_z;
  
  center_x = xDim / 2.0;
  center_y = yDim / 2.0;
  center_z = zDim / 2.0;

  float theta, phi;

  theta = RandomNG::randFloat(0.0, 2 * PI);
  phi = RandomNG::randFloat(0.0, 2 * PI);

  float x, y, z;

  x = nucleus_radius * cos(theta) * sin(phi) + center_x;
  y = nucleus_radius * sin(theta) * sin(phi) + center_y;
  z = nucleus_radius * cos(phi) + center_z;

  Point p;

  p.x = (int) x;
  p.y = (int) y;
  p.z = (int) z;

  return p;
}

Point LatticeEnvironment :: pickRandomPointOnCellMembrane()
{
  Point picked_point;

  picked_point.x = RandomNG::randInt(1, xDim - 1);
  picked_point.y = RandomNG::randInt(1, yDim - 1);
  picked_point.z = RandomNG::randInt(1, zDim - 1);

  int surface_area;

  surface_area
    = 2 * xDim * yDim
    + 2 * yDim * zDim
    + 2 * zDim * xDim;

  int random_for_face_axis;
  random_for_face_axis = RandomNG::randInt(1,surface_area);

  int random_for_face_sign;
  random_for_face_sign = RandomNG::randInt(1,2);

  if (random_for_face_axis <= 2 * xDim * yDim)
  {
    if (random_for_face_sign == 1)
    {
      picked_point.z = 0;
    }
    else
    {
      picked_point.z = zDim;
    }
  }
  else if (random_for_face_axis - 2 * xDim * yDim <= 2 * yDim * zDim)
  {
    if (random_for_face_sign == 1)
    {
      picked_point.x = 0;
    }
    else
    {
      picked_point.x = xDim;
    }
  }
  else
  {
    if (random_for_face_sign == 1)
    {
      picked_point.y = 0;
    }
    else
    {
      picked_point.y = yDim;
    }
  }

  return picked_point;
}

void LatticeEnvironment :: countAndReportAccessibleNucleus(){

  int exposed_nucleus_volume;
  int hidden_nucleus_volume = 0;

  exposed_nucleus_volume = countAndRecordAccessibleNucleus(NULL, NULL);

  Point * exposed_points = (Point *) malloc (exposed_nucleus_volume * sizeof (Point));

  countAndRecordAccessibleNucleus(exposed_points, NULL);

  while (exposed_nucleus_volume - hidden_nucleus_volume > NUCLEUS_VOLUME_TO_EXPOSE)
  {
    int i_point;

    i_point = RandomNG::randInt(0, exposed_nucleus_volume - 1);

    if (exposed_points[i_point].x >= 0)
    {
      arr[(int) exposed_points[i_point].x][(int) exposed_points[i_point].y][(int) exposed_points[i_point].z]->setEnvironment(INACCESSIBLE_ENVIRONMENT);
      exposed_points[i_point].x = -1;
      hidden_nucleus_volume ++;
    }
  }

  free (exposed_points);
  //temp
  //cout << "Exposed nucleus before:\t" << exposed_nucleus_volume << ", after: " 
  //     << ((NUCLEUS_VOLUME_TO_EXPOSE < exposed_nucleus_volume) ? NUCLEUS_VOLUME_TO_EXPOSE : exposed_nucleus_volume) << endl;
    
}

int LatticeEnvironment :: countAndRecordAccessibleNucleus(Point * accessible_nucleus_points, bool *** accessible_nucleus_lattice){
  int exposed_nucleus_volume = 0;
  int px, py, pz;

  for (px = 0; px < xDim; px ++)
    for (py = 0; py < yDim; py ++)
      for (pz = 0; pz < zDim; pz ++)
      {
	bool exposed = false;

	int env = arr[px][py][pz]->getEnvironment();

	if (env == CYTOSOL_ENVIRONMENT)
	{
	  int neighbour_x, neighbour_y, neighbour_z;

	  for (neighbour_x = px - 1; neighbour_x <= px + 1; neighbour_x ++)
	    for (neighbour_y = py - 1; neighbour_y <= py + 1; neighbour_y ++)
	      for (neighbour_z = pz - 1; neighbour_z <= pz + 1; neighbour_z ++)
	      {
		int neighbour_env;
		if (neighbour_x >= 0 && neighbour_x < xDim && neighbour_y >= 0 && neighbour_y < yDim && neighbour_z >= 0 && neighbour_z < zDim)
		{
		  neighbour_env = arr[neighbour_x][neighbour_y][neighbour_z]->getEnvironment();
		  if (neighbour_env == NUCLEUS_ENVIRONMENT)
		  {
		    exposed = true;
		  }
		}
	      }
	  if (exposed)
	  {
	    if (accessible_nucleus_points != NULL)
	    {
	      accessible_nucleus_points[exposed_nucleus_volume].x = px;
	      accessible_nucleus_points[exposed_nucleus_volume].y = py;
	      accessible_nucleus_points[exposed_nucleus_volume].z = pz;
	    }
	    if (accessible_nucleus_lattice != NULL)
	    {
	      accessible_nucleus_lattice[px][py][pz] = true;
	    }
	    exposed_nucleus_volume ++;
	  }
	}
      }

  return exposed_nucleus_volume;
}

#if 0
//////////////////////////////////////////////////////////////////////
//  Construct the environment for the Lattice space
//
void LatticeEnvironment :: buildLatticeEnvironment(){

	string line;
	string curLine;
	string delimiter = "\t";
	ifstream inFile (INPUT_FILE_ENV);

	while (!inFile.eof())
	{	
		getline(inFile,line,'\n');

		if(line.length() > 0){
			int x,y,z, env;
		
			int keyPos = 0;
			
			// Get X position
			keyPos = line.find(delimiter);	
			x = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get Y position
			keyPos = line.find(delimiter);	
			y = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get Z position
			keyPos = line.find(delimiter);	
			z = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get the environment
			env = atoi(line.substr(0).c_str());

			// Set the environment at position (x,y,z) to env
			arr[x][y][z]->setEnvironment(env);
			//printf("(%d,%d,%d) -> %d\n",x,y,z, arr[x][y][z]->getEnvironment());
			
		}
	}
	
	inFile.close();
}
#endif

//////////////////////////////////////////////////////////////////////
//  Construct an environment table that holds:
//  Environment type -> compatible environment (a list of vertices in lattice space)
//
void LatticeEnvironment :: buildEnvironmentTable(){

  //  map<moltype, vector <Molecule *> >::const_iterator itr = molTable.begin();
  //  for (itr = molTable.begin(); itr != molTable.end(); itr ++)
  //  {
  //  }

//  for (int i = 0; i < compatibleTable.size(); i ++)
//   map <moltype, vector <envtype> > compatibleTable;


	// Step through the array and for each Lattice Site
	for(int x = 0 ; x < xDim ; x ++){
		for(int y = 0 ; y < yDim ; y ++){
			for(int z = 0 ; z < zDim ; z ++){
			
				// The Lattice Site at location (x,y,z) of the Lattice
				LatticeSite * ls = arr[x][y][z];

				// Add its location to the environment table according to
				// its corresponding environment
				envTable[ls->getEnvironment()].push_back(ls->getLocation());

				// Build compatiblePointsTable
				map <moltype, vector <envtype> >::const_iterator itr;

				for (itr = compatibleTable.begin(); itr != compatibleTable.end(); itr ++){
				  //				  vector <envtype>::const_iterator env_itr;

				  if (isCompatibleEnvironment(itr->first, ls->getEnvironment())){
				    compatiblePointsTable[itr->first].push_back(ls->getLocation());
				  }
				}
			}
		}
	}


}



//////////////////////////////////////////////////////////////////////
//  Construct the list of environments that molecules are compatible with
//
void LatticeEnvironment::buildMoleculeCompatibleEnvTable(){

	string line;
	string curLine;
	string delimiter = "\t";
	ifstream inFile (INPUT_FILE_COMPATIBLE_ENV);

	while (!inFile.eof())
	{	
		getline(inFile,line,'\n');

		if(line.length() > 0){
			moltype mol;
			envtype env;
		
			int mol_num;
			int env_num;

			assert(sscanf (line.c_str(), "%d %d", &mol_num, &env_num) == 2);
			mol = (moltype) mol_num;
			env = (envtype) env_num;
			
			compatibleTable[mol].push_back(env);
		}
	}
	
	inFile.close();
	
}



//////////////////////////////////////////////////////////////////////
//  Set the LatticeSite object at the specified point in the Lattice Space
//
LatticeSite * LatticeEnvironment :: getLatticeSite(Point & p){

	if(withinCell(p)){
		// Lattice Site at Point p (Make sure its an INTEGER value)
		int x = roundToInt(p.x);
		int y = roundToInt(p.y);
		int z = roundToInt(p.z);

		// Return the pointer to the Latticesite at (x,y,z)
		return arr[x][y][z];
	}
	
	// LatticeSite is not within the cell, return NULL
	return (LatticeSite *)0;
}



//////////////////////////////////////////////////////////////////////
//  Get the Environment type at location p
//
int LatticeEnvironment :: getEnvironment(Point & p){

	LatticeSite * l = getLatticeSite(p);
	if((LatticeSite *)l != 0)
		return getLatticeSite(p)->getEnvironment();

	// location not within the cell
	return -1;
}



//////////////////////////////////////////////////////////////////////
//  Determine whether the specified molecule type is compatible in 
//  the specified environment
//
bool LatticeEnvironment :: isCompatibleEnvironment(moltype mol, envtype env){

	int size = compatibleTable[mol].size();

	// Step through the list of compatible environments for the molecule type 'mol'
	for(int envIndex = 0 ; envIndex < size ; envIndex ++){
	  // if the molecule type 'mol' is compatible with the environment 'env'
	  if(compatibleTable[mol][envIndex] == env)
	  {		
	    // Is a compatible environment
	    return true;
	  }
	}

	// Environment 'env' is not compatible for molecule type 'mol'
      	return false;
}



//////////////////////////////////////////////////////////////////////
//  Determine whether the provided location is within the defined dimensions
//  of the Lattice space
//
bool LatticeEnvironment :: withinCell(Point & p){

	int x = roundToInt(p.x);
	int y = roundToInt(p.y);
	int z = roundToInt(p.z);
	// Is specified Point within dimensions
	return 	withinCell(x,y,z);
}

bool LatticeEnvironment :: withinCell(int x, int y, int z){

	// Is specified Point within dimensions
        return  (x < xDim) && (x >= 0) &&
                        (y < yDim) && (y >= 0) &&
                        (z < zDim) && (z >= 0);

}

//////////////////////////////////////////////////////////////////////
//  Get a list of available locations for the specified environment
//
vector <Point *> * LatticeEnvironment :: getAvailableSites(moltype mol){

	// return a list of available locations with the environment 'env'
	return &(compatiblePointsTable[mol]);
}



//////////////////////////////////////////////////////////////////////
//  Display the Lattice Environments of each LatticeSite in the Lattice Space
//
void LatticeEnvironment :: displayLatticeEnv(){
	for(int x = 0 ; x < xDim ; x ++)
		for(int y = 0 ; y < yDim ; y ++)
			for(int z = 0 ; z < zDim ; z ++)
				// Output Location and its environment
				printf("(%d,%d,%d) -> %d\n"
					,x,y,z,(arr[x][y][z])->getEnvironment());
}


void LatticeEnvironment :: updateStateMatrix(){
	stateMatrix->add(updateMatrix);
}


void LatticeEnvironment :: populateStateMatrix(){
  //#define UNIFORM_CALCIUM
  //#define CORNER_CALCIUM_PULSE
	// this section is hardcoded; an external means of defining this would be
	// preferable in the longer term
#ifdef UNIFORM_CALCIUM
	for (int x = 0; x < xDim; x++) {
		for (int y = 0; y < yDim; y++) {
			for (int z = 0; z < zDim; z++) {
				int id = arr[x][y][z]->getID();
				Tuple * t = stateMatrix->getRow(id);
				t->tuple[0] = 1.0e-3; // This is 1.0 mM BPG
			}
		}
	}
#endif
#ifdef CORNER_CALCIUM_PULSE
	//	int id = arr[xDim/2][xDim/2][xDim/2]->getID();
	int id = arr[0][0][0]->getID();
	Tuple * t = stateMatrix->getRow(id);
	t->tuple[0] = 1.0e0;
#endif

	Tuple * curStat = stateMatrix->collapse();
        statisticsMatrix->appendRow(curStat);
}


void LatticeEnvironment :: displayStateMatrix(){
	for(int i = 0; i < numLatticeSites ; i ++){
		for(int j = 0 ; j < NUM_META ; j ++){
			printf("%f\t",stateMatrix->matrix[i]->tuple[j]);
		}
		printf("\n");
	}
}

// Build the table that specifies the diffusion rate of metabolite 'k' from env 'i' to env 'j'
// 0 <= Diffusion Rate <= 1 , fraction of diffusable metabolite independently
void LatticeEnvironment :: buildDiffusionRateTable(){
	diffusionRateTable = new Tuple ** [numEnvironments];
	for(int i = 0; i < numEnvironments ; i ++){
		diffusionRateTable[i] = new Tuple * [numEnvironments];
	}

        string line;
        string curLine;
        string delimiter = "\t";
        ifstream inFile (INPUT_FILE_DIFFUSION_RATES);

        while (!inFile.eof())
        {
                getline(inFile,line,'\n');

                if(line.length() > 0){
                        envtype from, to;

                        int keyPos = 0;

                        // Get X position
                        keyPos = line.find(delimiter);
                        from = atoi(line.substr(0,keyPos).c_str());
                        line = line.substr(keyPos+delimiter.length());

                        // Get Y position
                        keyPos = line.find(delimiter);
                        to = atoi(line.substr(0,keyPos).c_str());
                        line = line.substr(keyPos+delimiter.length());

			diffusionRateTable[from][to] = new Tuple(NUM_META);

			for(int i = 0 ; i < NUM_META ; i ++){
                        	keyPos = line.find(delimiter);
	                        diffusionRateTable[from][to]->tuple[i] = atof(line.substr(0,keyPos).c_str());
        	                line = line.substr(keyPos+delimiter.length());
			}
                }
        }

        inFile.close();
/*
	for(int i = 0 ; i < numEnvironments ; i ++){
		diffusionRateTable[i] = new Tuple * [numEnvironments];
		for(int j = 0 ; j < numEnvironments ; j ++){
			diffusionRateTable[i][j] = new Tuple(NUM_META);
			for(int k = 0 ; k < NUM_META ; k ++){
				diffusionRateTable[i][j]->tuple[k] = 0.8 / (k*0.25+1);
			}
		}
	}
*/
}

void LatticeEnvironment :: buildReactionRateTable(){
	reactionRateTable = new Tuple * [numEnvironments];
	for(int i = 0 ; i < numEnvironments ; i ++){
		reactionRateTable[i] = new Tuple(NUM_MOL_ENZYME);
		for(int j = 0 ; j < NUM_MOL_ENZYME ; j ++){
			reactionRateTable[i]->tuple[j] = 0.5;
		}
	}
}

struct diffuseInto {
	float			multiplier;
	vector<int>  		ids;
	struct diffuseInto *	next;
};

struct diffuseInto *** diffusionPatterns;

void LatticeEnvironment :: buildDiffusionPatterns(){
	int i, k, m, met;

	diffusionPatterns = new diffuseInto **[NUM_META];

	for (i = 0; i < NUM_META; i++) {
		diffusionPatterns[i] = new diffuseInto *[xDim*yDim*zDim];
		for (int k = 0; k < xDim*yDim*zDim; k++)
			diffusionPatterns[i][k] = NULL;
	}

	//#define ONE_AXIS_DIFFUSION

	float ratio[3];
	ratio[0] = 1;//1.0;
#ifdef ONE_AXIS_DIFFUSION
	ratio[1] = 0;
	ratio[2] = 0;
#else
	ratio[1] = 0.707107;//sqrt(2.0);
	ratio[2] = 0.57735;//sqrt(3.0);
#endif /* ONE_AXIS_DIFFUSION */

	//	ratio[0] = .1;
	//	ratio[1] = 0;
	//	ratio[2] = 0;

	float ratio_sum;

	ratio_sum = 6 * ratio[0] + 12 * ratio[1] + 8 * ratio[2];

	ratio[0] *= 6 / ratio_sum;
	ratio[1] *= 12 / ratio_sum;
	ratio[2] *= 8 / ratio_sum;

	// Number of Dimensions
	int numDegrees = 3;

	// Step through each lattice site in the system
	for(i = 0 ; i < xDim ; i ++){
	for(k = 0 ; k < yDim ; k ++){
	for(m = 0 ; m < zDim ; m ++){

		LatticeSite * origin = arr[i][k][m];

		// ID of the source of diffusion
		int origin_id = origin->getID();

		// Get the source environment
		envtype fromEnv = origin->getEnvironment();

		float remove[NUM_META];
		for (met = 0; met < NUM_META; met++) remove[met] = 0.0;

		/////////////////////////////////////////////////////////////////////
		// Calculate diffusion rates to neighbors with each degree
		for(int degree = 0 ; degree < numDegrees ; degree++){

		// Get the neighbors of the source lattice
		vector<LatticeSite *> * neighborLattices = &(origin->neighbors[degree]);
		int numLattices = neighborLattices->size();

		for (int num = 0; num < numLattices; num++) {
			// Get the neighboring lattice site
			LatticeSite * neighbor = (*neighborLattices)[num];
			int neighbor_id = neighbor->getID();
			// determine the neighboring lattice's environment
			envtype toEnv = neighbor->getEnvironment();
			// determine its corresponding diffusion rate (origin->destination env)
			Tuple * neighborRate = diffusionRateTable[fromEnv][toEnv];

			for (met = 0; met < NUM_META; met++) {
				float multiplier = neighborRate->tuple[met];
				if (multiplier <= 0.0) continue;

				multiplier *= ratio[degree] * TIMESCALE / (SPACESCALE * SPACESCALE);
				//				cout << multiplier << " = " << neighborRate->tuple[met] << " * " << ratio[degree] << " * " << TIMESCALE << " * " << SPACESCALE << "^2" << endl;
				assert (multiplier <= 0.1);

				remove[met] -= multiplier;

				diffuseInto * dif = diffusionPatterns[met][origin_id];

				if (dif == NULL) {
					dif = new diffuseInto;
					dif->multiplier = multiplier;
					dif->ids.push_back(neighbor_id);
					dif->next = NULL;

					diffusionPatterns[met][origin_id] = dif;
				}
				else while (1) {
	                                if (dif->multiplier == multiplier) {
						dif->ids.push_back(neighbor_id);
						break;
					}

					if (!dif->next) {
						diffuseInto * dif2 = new diffuseInto;
						dif2->multiplier = multiplier;
						dif2->ids.push_back(neighbor_id);
						dif2->next = NULL;

						dif->next = dif2;

						break;
					}

					dif = dif->next;
				}

			}
		}
		}

		for (met = 0; met < NUM_META; met++) {
			diffuseInto * dif = new diffuseInto;
			dif->multiplier = remove[met];
			dif->ids.push_back(origin_id);
			dif->next = diffusionPatterns[met][origin_id];

			diffusionPatterns[met][origin_id] = dif;
		}
	}
	}
	}
}

void LatticeEnvironment :: diffuse() {
	int maxIdx = xDim * yDim * zDim;

#if 0	
	{
	  static int timestep = 0;
#if 0
	  int x, y;
	  for (x = 0; x < xDim; x ++)
	  {
	    cout << timestep <<",,";
	    for (y = 0; y < yDim; y ++)
	    {
	      int id = arr[x][y][xDim / 2]->getID();
	      Tuple * row = stateMatrix->getRow(id);
	      cout << row->tuple[0] << ",";
	    }
	    cout << endl;
	  }
#endif

	  int id = arr[xDim / 2][yDim / 2][xDim / 2]->getID();
	  Tuple * row = stateMatrix->getRow(id);
	  cout << timestep << ",\t" << row->tuple[0] << ",\t" << row->tuple[1] << ",\t" << row->tuple[2] << ",\t" << row->tuple[3] << ",\t" << row->tuple[4] << endl;

	  timestep ++;
	}
#endif

	for (int idx = 0; idx < maxIdx; idx++) {
		Tuple * t = stateMatrix->getRow(idx);

		for (int met = 0; met < NUM_META; met++) {
			float conc = t->tuple[met];
			if (conc <= 0.0) continue;

			float removed = 0.0;

			diffuseInto * dif = diffusionPatterns[met][idx];

			while (dif) {
				int num = dif->ids.size();

				for (int i = 0; i < num; i++) {
					int neighbor_id = dif->ids.at(i);
					Tuple * ut = updateMatrix->getRow(neighbor_id);
					ut->tuple[met] += conc * dif->multiplier;
					removed += conc * dif->multiplier;
				}

				dif = dif->next;
			}

			t->tuple[met] -= removed;
		}
	}

	// Update the state of the Cell (Concentrations)
	stateMatrix->addAndClear(updateMatrix);
}
#if 0
void LatticeEnvironment :: buildDiffusionPatterns(){}

void LatticeEnvironment :: diffuse(){

	float ratio[3];
	ratio[0] = 1;//1.0;
	ratio[1] = 0.707107;//sqrt(2.0);
	ratio[2] = 0.57735;//sqrt(3.0);

	// Step through each lattice site in the system
	for(int i = 0 ; i < xDim ; i ++){
	for(int j = 0 ; j < yDim ; j ++){
	for(int k = 0 ; k < zDim ; k ++){

		LatticeSite * origin = arr[i][j][k];

		// Number of Dimensions
		int numDegrees = 3;

		// ID of the source of diffusion
		int origin_id = origin->getID();

		// Get the source environment
		envtype fromEnv = origin->getEnvironment();

		// Table to store diffusion rates to neighboring lattices
		Tuple ** neighborRates;
		neighborRates = new Tuple *[numDegrees];

		// table of net diffusion matrices, one for each degree
		Matrix ** diffusionRates;
		diffusionRates = new Matrix *[numDegrees];

		// update order for each degree of neighbors
		vector<int> updateOrder[numDegrees];

		// Tuple representing total diffusion rate of the metabolites
		Tuple * totalRates = new Tuple(NUM_META);

		/////////////////////////////////////////////////////////////////////
		// Calculate diffusion rates to neighbors with each degree
		for(int degree = 0 ; degree < numDegrees ; degree++){

		// Get the neighbors of the source lattice
		vector<LatticeSite *> * neighborLattices = &(origin->neighbors[degree]);
		int numLattices = neighborLattices->size();
		// create a diffusion matrix where # rows = # neighbors
		diffusionRates[degree] = new Matrix(numLattices,NUM_META);

		// Iterate each of the neighbors
		for(int num = 0 ; num < numLattices ; num ++){
			// Get the neighboring lattice site
			LatticeSite * neighbor = (*neighborLattices)[num];
			int neighbor_id = neighbor->getID();
			// add to the update order list
			updateOrder[degree].push_back(neighbor_id);
			// determine the neighboring lattice's environment
			envtype toEnv = neighbor->getEnvironment();
			// determine its corresponding diffusion rate (origin->destination env)
			Tuple * neighborRate = diffusionRateTable[fromEnv][toEnv];
			// insert the diffusion rate to the diffusion matrix
			diffusionRates[degree]->copyValueToRow(num,neighborRate);
		}

		// Calculate the total relative diffusion rates to all neighbors
		neighborRates[degree] = diffusionRates[degree]->collapse();
		// diffusion rate varies according to its degree
		neighborRates[degree]->scalarize(ratio[degree]);

		// add diffusion rate of the specified degree to the overall diffusion rate
		totalRates->add(neighborRates[degree]);
		delete neighborRates[degree];
		}
		/////////////////////////////////////////////////////////////////////
		
		/////////////////////////////////////////////////////////////////////
		// Caluculate total diffusable ratio
		Tuple * updateOriginConc = stateMatrix->getRow(origin_id);
		Tuple * netRates = new Tuple(totalRates);
		netRates->scalarize(1.0/26.0);
		/////////////////////////////////////////////////////////////////////

		/////////////////////////////////////////////////////////////////////
		// Update Neighbor concentrations on the update matrix
		for(int degree = 0 ; degree < numDegrees ; degree ++){
			int num = diffusionRates[degree]->getNumRows();
			for(int drate = 0 ; drate < num ; drate ++){
				Tuple * t = new Tuple(diffusionRates[degree]->getRow(drate));
				t->scalarize(ratio[degree]/26.0);
				t->divByElement(totalRates);
				t->multByElement(netRates);
				t->multByElement(updateOriginConc);

				updateMatrix->addValueToRow(updateOrder[degree][drate],t);
				delete t;
			}
			delete diffusionRates[degree];
		}
		/////////////////////////////////////////////////////////////////////
		// Update the concentration of the original site
		netRates->negate();

		netRates->multByElement(updateOriginConc);
                updateMatrix->addValueToRow(origin_id,netRates);

		// clean up
		delete netRates;
		delete totalRates;
		delete [] diffusionRates;
		delete [] neighborRates;		

	}
	}
	}

	// Update the state of the Cell (Concentrations)
	stateMatrix->addAndClear(updateMatrix);
}
#endif

void LatticeEnvironment :: pumpInMetabolite(float amount_to_pump_in){
	int id = arr[0][0][0]->getID();
	Tuple * row = stateMatrix->getRow(id);
	row->tuple[META_A] += amount_to_pump_in;
}

void LatticeEnvironment :: pumpOutMetabolite(float amount_to_pump_out){
	int id = arr[xDim - 1][yDim - 1][zDim - 1]->getID();
	Tuple * row = stateMatrix->getRow(id);
	row->tuple[4] -= amount_to_pump_out;
	if (row->tuple[META_E] < 0)
	{
	  row->tuple[META_E] = 0;
	}
}


void LatticeEnvironment :: reactWithMetabolite(Point & point, int reactant, int product, int enzType, float forward_rate, float reverse_rate, float equilibrium_constant, float maximum_rate){
	LatticeSite * l = getLatticeSite(point);
	int id = l->getID();

	float amount;

	float Vm, a, p, Keq, Ka, Kp;

	Vm = maximum_rate * TIMESCALE;
	a = stateMatrix->get(id,reactant);
	p = stateMatrix->get(id,product);
	Keq = equilibrium_constant;
	Ka = forward_rate;
	Kp = reverse_rate;

	const static float MINIMUM_CONCENTRATION = 1.0e-10;

	//		cout << "Vm = " << Vm << ", a = " << a << ", p = " << p << ", Keq = " << Keq << ", Ka = " << Ka << ", Kp = " << Kp << endl;

	if (a >= MINIMUM_CONCENTRATION)
	{
	  if (p >= MINIMUM_CONCENTRATION)
	  {
	    amount = Vm * ((a/Ka)*(1-((p/a)/Keq))) / (1 + (a/Ka) + (p/Kp));
	  }
	  else
	  {
	    amount = Vm * (a/Ka) / (1 + (a/Ka));
	  }
	  //		  forward_amount = Vm * ((a/Ka)*(1-((MAR)/Keq))) / (1 + (a/Ka) + (p/Kp));
	}
	else
	{
	  amount = 0.0;
	}

	double cube_volume;

	cube_volume = SPACESCALE * SPACESCALE * SPACESCALE; /* volume in cubic metres */;
	cube_volume *= 1000; /* volume in litres */
#if 0
	cout << "conc a = " << a << ", conc p = " << p << ". ";

	cout << "a (" << a << ") / Ka (" << Ka << ") = " << a / Ka << ". ";
	cout << "p (" << p << ") / Kp (" << Kp << ") = " << p / Kp << ". ";
	cout << "p/a (" << p/a << ") / Keq (" << Keq << ") = " << (p/a) / Keq << ". ";
	cout << "Vm = " << Vm << ". ";

	cout << "Converting " << amount << " moles: ";
#endif
	amount /= cube_volume; /* divide by VOLUME_OF_CUBE to get concentration change */;

	//		cout << amount << " mol/L concentration change in the " << SPACESCALE << " on-a-side cube." << endl;
		
	if (amount > a)
	{
	  amount = a;
	}
	else if (-amount > p)
	{
	  amount = -p;
	}

	/* This is equation (A1) from page 5326 of Eur. J. Biochem. 267, (2000).
	 * It is from the article "Can yeast glycolysis be understood in terms of in vitro
	 * kinetics of the constituent enzymes? Testing biochemistry."
	 */

	if(amount == 0)	return;

	Tuple * t = new Tuple(NUM_META);
	t->tuple[reactant] = -1 * amount;
	t->tuple[product] = amount;
	stateMatrix->addValueToRow(id,t);

	delete t;
}

void LatticeEnvironment :: displayState(){
	stateMatrix->display();
}

Matrix * LatticeEnvironment :: getStateMatrix(){
	return 	stateMatrix;
}

Matrix * LatticeEnvironment :: getStatisticsMatrix(){
	return 	statisticsMatrix;
}

LatticeMacroFeature * LatticeEnvironment :: getLatticeMacroFeatures(){
	return 	latticeMacroFeatures;
}

Tuple * LatticeEnvironment :: getCurStat(){
	int num = statisticsMatrix->getNumRows();
	if(num > 0)
		return statisticsMatrix->matrix[num-1];
	return (Tuple *)0;
}

int LatticeEnvironment :: firstEnv(moltype mol) {
  return compatibleTable[mol][0];
}
