// start break
// beginning
//package org.opensourcephysics.sip.ch08.md;
import java.awt.*;
import org.opensourcephysics.display.*;
import org.opensourcephysics.frames.*;
import org.opensourcephysics.numerics.*;

/**
 * LJParticles simulates a two-dimensional system of particles which interact via a Lennard-Jones potential.
 *
 * Through user inputs from LJParticlesApp, these particles can be made to be a granular media.
 *
 * Daniel Costantino, based on LJParticles by:
 * @author Jan Tobochnik, Wolfgang Christian, Harvey Gould
 * @version 1.0 revised 03/28/05
 */

public class LJParticles implements Drawable, ODE {

   public double state[];	//state array for the particles
   public double ax[], ay[];	//x and y-componenets of acceleration
   //public int graphArray[];
   public int N, nx, ny;	// number of particles, number per row, number per column (nx and ny were used for initial
				// conditions in the original LJParticles but are no longer used)
   public int neighbors[][];	// neighbors keeps track of the near neighbors of each particle
   public int numberNeighbors[];// numberNeighbors keeps track of how many neighbors each parictle has
   public double Lx, Ly;	// the width and height of the particles container
   public double rho = N/(Lx*Ly);	// the particle density is the number of particles / area
   public double initialKineticEnergy;	// the initial kinetic energy per particle
   public double rCutOff, rCutOffSq;	// the cut-off distance beyond which particles will not feel any force between them
					// to eliminate attractive forces which are unphysical for granular matter
   public double gamma;		// coefficient of a damping factor;
   public int steps = 0;	// keeps track of how many time steps have passed to correctly average energies
   public int stepsBetween = 0;	// a user input to manually determine how often the neighbors list should be updated
   public double dt = 0.01;	// how long each time step is
   public double inputDt;	// the time step the user chooses initially (for some methods in keeping the temperature
				// under control, this is used as a comparison
   public double t;		// the current time
   public double computeTime;	// keeps track of when neighbors list should be calculated (using stepsBetween * dt)
   public double computeTime2;	// keeps track of when neighbors list should be updated (using calculation of how far particles
				// should have been able to move
   //public double dampingV2, damper;	// inputs for a drastic damping measure which isn't being used since it didn't work right
   public double totalPotentialEnergyAccumulator; // the total potential energy of the system
   public double totalKineticEnergyAccumulator, totalKineticEnergySquaredAccumulator; // the total kinetic energy and KE^2
					// of the system
   public double virialAccumulator;	// the virial of the system (has no meaning at this point, but I didn't want to remove
					// the book authors' work in case it becomes useful later)
   //public String initialConfiguration;	// if set-ups other than random are allowed
   public double radius = 0.5; // radius of particles on screen
   public double phi = 0;	// the angle giving the particles current direction of motion (for some methods of controlling
				// temeperature)
   public double desiredKE;	// the maximum KE per particle which is added via random kicks
   public boolean highKE = false;	// whether or not a particle has a kinetic energy high above the initial value (for
				// controlling the temperature
   public double lastKineticEnergy = 0;	// keeps track of what the kinetic energy was for the system on the last pass
   public boolean smallR[];	// keeps track of if there are any particles very close together
   public boolean allBigR = true;	// keeps track of if there are any particles very close together
   public boolean lastSmall[];	// keeps track of if there are any particles very close together
   public boolean giveRandom = false;	// whether or not random kicks should be given to the particles
   public boolean nuclear = false;	// whether or not a "nuclear" option in reducing KE should be used
   public double holeDiam;	// the size of a hole in a silo simulation
   public int nPart;		// how many particles there are at any given time.
   public double g = 9.81;	// acceleration due to gravity
   public int consider[];	// whether or not a particluar particle should be looked at (no longer used, but could be
				// useful in the future)
   public boolean isSilo = false;	// whether or not this is a silo simulation
   //public double graphTime;
   public boolean bootstrap = false;	// whether or not to use a drastic measure of temperature reduction
   public boolean theMatrix = false;	// whether or not to allow time to slow down if particles get too close together

   Verlet verlet_solver = new Verlet(this);	// opening a Verlet solver for the system

//initialize() initializes the system---------------------------------------------------------------------------------------
   public void initialize() {
      //N = nx*ny;
      nPart = N; 		// the initial nPart should be N
      t = 0;			// initial time is 0
      rho = N/(Lx*Ly);		// density (which is constant for the current simulations)
      resetAverages();		// set all energies to zero
      state = new double[1+4*N];// each of N particles gets 4 pieces of state to track it; for the i-th particle:
				// state[4 * i]  = x
				// state[4 * i + 1] = vx
				// state[4 * i + 2] = y
				// state[4 * i + 3] = vy
      consider = new int[N];	// an array to see whether or not particles should be looked at

      lastSmall = new boolean[N];
      smallR = new boolean[N];

	int k = 0;		// counter for use in initializing consider:
		while(k < N) {
			consider[k] = 1;
			lastSmall[k] = false;
			smallR[k] = false;
			k = k + 1;
		}
      ax = new double[N];	// acceleration in x-direction
      ay = new double[N];	// acceleration in y-direction
    					  //if(initialConfiguration.equals("crystal")) {
      						  // setTriangularLattice();
     					 //} else if(initialConfiguration.equals("rectangular")) {
        					// setRectangularLattice();
     					 //} else {

         setRandomPositions();	// it really only makes sense to allow for random positions, but I left the others in incase
				// I change my mind
      					//}
      setVelocities();		// get the initial particle velocities

      neighbors = new int[N][N];// the neighbors list should allow for each particle to neighbor every other (as an extreme case)
				// the second size really only needs to be N - 1
      numberNeighbors = new int[N];//again this size should really only need to be N - 1
      rCutOffSq = rCutOff * rCutOff;// finding the square of the cut-off distance  (cut-off distance is a user input)

      computeTime2 =  (2.7 - rCutOff) / (2 * Math.sqrt(2*initialKineticEnergy));
				// theoretical calculation of how often user list should be calculated
	computeTime = stepsBetween * dt;	// allows for user input to decide this
	//System.out.println("compute Time = "+computeTime+"   compute Time2 =    "+computeTime2);

	inputDt = dt;
      findNeighbors();		// see which particles are near which other particles
      computeAcceleration();	// find the initial accleration

      verlet_solver.setStepSize(dt);	// prepare the Verlet algorithm
   }

//setRandomPositions() places the particles at random positions in the container, but no closer than the value
//rMinimumSquared----------------------------------------------------------------------------------------------------------
   public void setRandomPositions() {
      double rMinimumSquared = Math.pow(2.0, 1.0/3.0);	// to keep particles from getting right on top of each other
      boolean overlap;					// same reasoning
      for(int i = 0; i<N; ++i) {
         do {
            overlap = false;
            state[4*i] = Lx*Math.random();   // get x and
            state[4*i+2] = Ly*Math.random(); // y positions
            int j = 0;
            while((j<i) && !overlap) {
               double dx = state[4*i]-state[4*j];// look at other particles
               double dy = state[4*i+2]-state[4*j+2];
               if(dx*dx+dy*dy < rMinimumSquared) {
                  overlap = true;	// if they're too close, then note that
               }
               j++;
            }
         } while(overlap);	// this should hang the system if particles are too close
      }
   }

// The other position functions are just commented out since I don't think I'll be using them for this project
   // end break
   // start break
   // setRectangularLattice
/*   public void setRectangularLattice() {     // place particles on a rectangular lattice
      double dx = Lx/nx; // distance between columns
      double dy = Ly/ny; // distance between rows
      int i = 0;
      int iy = 0;
      while(i<N) {
         for(int ix = 0; ix<nx; ++ix) { // loop through particles in a row
            if(i<N) {
               state[4*i] = dx*(ix+0.5);
               state[4*i+2] = dy*(iy+0.5);
               i++;
            }
         }
         iy++;
      }
   }

   // end break
   // start break
   // setTriangularLattice
   public void setTriangularLattice() { // place particles on triangular lattice
      double dnx = Math.sqrt(N);
      int ns = (int) dnx;
      if(dnx-ns>0.001) {
         ns++;
      }
      double dx = Lx/ns;
      double dy = Ly/ns;
      int i = 0;
      int iy = 0;
      while(i<N) {
         for(int ix = 0; ix<ns; ++ix) {
            if(i<N) {
               state[4*i+2] = dy*(iy+0.5);
               if(iy%2==0) {
                  state[4*i] = dx*(ix+0.25);
               } else {
                  state[4*i] = dx*(ix+0.75);
               }
               i++;
            }
         }
         iy++;
      }
   }*/

//setVelocities() gets the velocity of each particle, and then makes sure the center of mass's momentum is zero and
//rescales the total kinetic energy to match the user input----------------------------------------------------------------
   public void setVelocities() {
      double vxSum = 0.0;	// keep track of the total x and y-velocity of all particles
      double vySum = 0.0;
      for(int i = 0; i<N; ++i) {           // assign random initial velocities to each particle
         state[4*i+1] = Math.random()-0.5; // x and
         state[4*i+3] = Math.random()-0.5; // y components
         vxSum += state[4*i+1];		   // keep track of total x and y-velocity of all particles
         vySum += state[4*i+3];
      }

      double vxcm = vxSum/N; // find the center of mass momentum (velocity)
      double vycm = vySum/N;
      for(int i = 0; i<N; ++i) {// subtract this velocity from each particle's velocity
         state[4*i+1] -= vxcm;
         state[4*i+3] -= vycm;
      }

      double v2sum = 0; 
      for(int i = 0; i<N; ++i) {// find the total kinetic energy of the system by finding sum(v^2)
         v2sum += state[4*i+1]*state[4*i+1]+state[4*i+3]*state[4*i+3];
      }
      double kineticEnergyPerParticle = 0.5*v2sum/N;// then multiply v^2 by (m / 2) but m = 1
      double rescale = Math.sqrt(initialKineticEnergy/kineticEnergyPerParticle); // find the ratio between the kinetic
					// energy which whe have and what we want (from input)
      for(int i = 0; i<N; ++i) {	// force the velocities to get us the desired KE per particle
         state[4*i+1] *= rescale;
         state[4*i+3] *= rescale;
      }
   }

//From the total kinetic energy over all time, find the average kinetic energy (temperature) for 
   public double getMeanTemperature() {
      return totalKineticEnergyAccumulator/(N*steps);
   }

//Output the average energy of the system (averaged over time)
   public double getMeanEnergy() {
      return totalKineticEnergyAccumulator/steps+totalPotentialEnergyAccumulator/steps;
   }

/*
Pressure and Heat Capacity don't make any sense for my simulation, but I'm leaving them here in case I decide to play with
them later
   public double getMeanPressure() {
      double meanVirial;
      meanVirial = virialAccumulator/steps;
      return 1.0+0.5*meanVirial/(N*getMeanTemperature()); // quantity PA/NkT
   }

   public double getHeatCapacity() {
      double meanTemperature = getMeanTemperature();
      double meanTemperatureSquared = totalKineticEnergySquaredAccumulator/steps;
      // heat capacity related to fluctuations of kinetic energy
      double sigma2 = meanTemperatureSquared-meanTemperature*meanTemperature;
      double denom = sigma2/(N*meanTemperature*meanTemperature)-1.0;
      return N/denom;
   }
*/

//resetAVerages() sets all the energies equal to zero----------------------------------------------------------------------
   public void resetAverages() {
      steps = 0;
      virialAccumulator = 0;
      totalPotentialEnergyAccumulator = 0;
      totalKineticEnergyAccumulator = 0;
      totalKineticEnergySquaredAccumulator = 0;
   }

//computeAccleration() uses the neighbor list to look at nearby particles and find each particles acceleration
//-------------------------------------------------------------------------------------------------------------------------
   public void computeAcceleration() {
      int j = 0;
      int k = 0;
      for(int i = 0; i<N; i++) {// each time the acceleration of each particle should be set to zero
         ax[i] = 0;
         ay[i] = 0;
      }

      int randParticle = (int)(N * Math.random());	//choose a particle to kick at random

	double dx = 0;	// separation between i-th and j-th particles in x
	double dy = 0;	// and y
	double r2 = 0;	// separation^2
	double r = 0;	// sparation
	double fx = 0;	// force in x-
	double fy = 0;	// and y-directions
	double friction = 0;	// damping force
	double vx = 0;	// velocity in x
	double vy = 0;	// and y-directions
	double oneOverR2 = 0;	// 1/separation^2
	double oneOverR6 = 0;	// 1/separation^6
	double oneOverRCO2 = 1.0 / rCutOffSq;  //same but for cut-off distance
	double oneOverRCO6 = oneOverRCO2 * oneOverRCO2 * oneOverRCO2;
	double v2 = 0;	// velocity squared
	//double dampingV2 = 0;	// a damping force proportional to velocity which doesn't seem to work
	double fOverR = 0;	//force divided by R
	double fatRC = 0;	// force the particles would feel at the cut-off separation
	double bootStrapV = 0;	// a damping technique which doesn't seem to work

      for(int i = 0; i<=N-1; i++) {

         while(k < numberNeighbors[i])  {
			// cycle through all of i's neighbors (if it has any)
	    j = neighbors[i][k];
            dx = pbcSeparation(state[4*i]-state[4*j], Lx);	// find their separation in the x-
            dy = pbcSeparation(state[4*i+2]-state[4*j+2], Ly);	// and y-directions using periodic boundary conditions
            r2 = dx*dx+dy*dy;					// find the square of their separation
	    r = Math.sqrt(r2);					// and that separation
	    oneOverR2 = 1.0/r2;					// find 1/r^2
	    oneOverR6 = oneOverR2*oneOverR2*oneOverR2;		// and 1/r^6
	    v2 = (state[(4 * i) + 1] * state[(4 * i) + 1]) + (state[(4 * i) + 3] * state[(4 * i) + 3]);
								// find the i-th particles velocity^2
	    //dampingV2 = damper *  v2 / r;			// find this force to damp out high velocities, again
								// it doesn't seem to work

	    if( (r2 > rCutOffSq) || (nuclear && ((Math.abs(dx) < 0.15) && (Math.abs(dy) < 0.15)) ) ) {
		// for particles outside of the cut-off distance, there should be no force felt
		// or within too close a range if the "nuclear" option is selected
			fx = 0;
			fy = 0;
			friction = 0;
	    } else {
		fOverR = 48.0*oneOverR6*(oneOverR6-0.5)*oneOverR2;	// calculate the force felt between these particles
									// as described in the theory section of the guide
									// using epsilon and sigma as one
		fatRC = 48 * oneOverRCO6 * (oneOverRCO6 - 1/2) * oneOverRCO2;
		fx = (fOverR - fatRC) * dx;
		fy = (fOverR - fatRC) * dy;
		vx = state[(4 * i) + 1] - state[(4 * j) + 1];		// finding the relative velocities between the two
		vy = state[(4 * i) + 3] - state[(4 * j) + 3];		// particles i and j
		friction = -gamma * (vx * dx + vy * dy) / r2;		// finding the frictional force between them
		totalPotentialEnergyAccumulator += 4.0*((oneOverR6*oneOverR6-oneOverR6) - (oneOverRCO6 * oneOverRCO6 - oneOverRCO6));
					// to avoid discontinuity in energy calculate the potential energy using method
					// descibed in guide
	    }

	    if(theMatrix && (Math.abs(dx) < 0.15 && Math.abs(dy) < 0.15) && (dt > (inputDt / 50)) ) {
		dt = dt / 100;		// if the particles are too close together, reduce the time step (again
		smallR[i] = true;		// this doesn't seem to be helping much, but I think it's from where I
		allBigR = false;	// initialize allBigR and smallR)
	        verlet_solver.setStepSize(dt);	// prepare the Verlet algorithm
		//System.out.println("i "+i+" j "+j+" dx "+dx+" dy "+dy+" dt "+dt);
		//System.out.println("smallR[i] "+smallR[i]);
	    }

	    ax[i] += fx + (friction * dx)/*- dampingV2 * state[4 * i] * (-vXdir)*/;
            ay[i] += fy + (friction * dy)/*- dampingV2 * state[(4 * i) + 2] * (-vYdir)*/;
            ax[j] -= (fx + (friction * dx));
	    ay[j] -= (fy + (friction * dy));		// Calculate accleration using the mass as one and Newton's third
							// law

	    virialAccumulator += dx*fx+dy*fy;		// find the virial, which doesn't get used for anything at this point
	    k = k + 1;
         } // ends while(k < numberNeighbors[i])

		if((lastSmall[i] && !smallR[i]) && (dt != inputDt)) {
			//System.out.println("Here smallR[i] = "+smallR[i]);
			dt = inputDt;		// fixing the dt size if it is too small and all the particles are far away
						// from each other
		        verlet_solver.setStepSize(dt);	// prepare the Verlet algorithm
		}

		lastSmall[i] = smallR[i]; //if this one is small, then I next time I should fix the size
		smallR[i] = false;	// should assume that no r's are small next time

	if( giveRandom && i == randParticle && desiredKE != 0) {
		//Give particles a randomly sized and directed kick to give the added velocities a Maxwell Boltzman
		//distribution
		double rand = Math.random();
		double a = -Math.log(rand);
		double theta = 2.0 * Math.PI * Math.random();
		state[4 * i + 1] = Math.sqrt(2.0 * desiredKE * a) * Math.cos(theta);
		state[4 * i + 3] = Math.sqrt(2.0 * desiredKE * a) * Math.sin(theta);
	}
      } // ends for(i < N - 1)

   }

//pbcSeparation(dx, L) finds the separation between two particles using periodic boundary conditions
//-------------------------------------------------------------------------------------------------------------------------
   private double pbcSeparation(double ds, double L) {
      if(ds > 0) {	// with periodic boundary conditions, particles can be further apart than half the lengh of the
			// container
         while(ds > 0.5*L) {
            ds -= L;
         }
      } else {
         while(ds<-0.5*L) {
            ds += L;
         }
      }
      return ds;
   }

//pbcPosition(s, L) finds a particle's position using periodic boundary conditions
//------------------------------------------------------------------------------------------------------------------------
   private double pbcPosition(double s, double L) {
      if(s>0) {		// each particles postion must be in 0 < s < L
         while(s>L) {
            s -= L;
         }
      } else {
         while(s<0) {
            s += L;
         }
      }
      return s;
   }

//getRate(state, rate) uses the current state array and the Verlet solver to find the call computeAcceleration and
// findNeighbors------------------------------------------------------------------------------------------------------------
   public void getRate(double[] state, double[] rate) {
      // computeAcceleration() is called every other time getRate begins since the new velocities are calculated using
      // using the lastest two accelerations and in each step of the Verlet algorithm, getRate is called twice

      if( (t > computeTime) || (t > computeTime2) || (highKE && !isSilo)) {
		// if it's been a wahile, recalculate the neighborlist
	findNeighbors();

	computeTime2 = t + (2.7 - rCutOff) / (2 * Math.sqrt(2 * initialKineticEnergy));
	computeTime = t + stepsBetween * dt;	// when the neighbor list should be recalculated
	/*if(computeTime > computeTime2) {
		// output when things will be calculated again
		System.out.println("time = "+t+" computeTime = "+computeTime+" computeTime2 = "+computeTime2);
	} else {
		System.out.println("time = "+t+"  computeTime = "+computeTime);
	}*/
	highKE = false;
      }


      if(verlet_solver.getRateCounter()==1) {
         computeAcceleration();
	 /*if(dt != inputDt) {
		System.out.println("dt = "+dt);
	 }*/
	 //smallR = false;
	 allBigR = true; 

      }

		
      for(int i = 0; i<N; i++) {
         rate[4*i] = state[4*i+1];   // rates for positions are velocities
         rate[4*i+2] = state[4*i+3];
         rate[4*i+1] = ax[i];        // rate for velocity is acceleration
         rate[4*i+3] = ay[i];
      }
      //if(dt == inputDt) {
	      rate[4*N] = 1; // since dt/dt = 1
      //} else {
	//      rate[4 * N] = inputDt / dt;
      //}//I'm not sure if this would work or not, and since everything seems fine at this point, I'm going with if it's not
	 // broke don't fix it
   }

//getstate() calculates the state of the system via the Verlet method------------------------------------------------------
   public double[] getState() {
      return state;
   }

//at each step the histograms for x and y-velocities and speed are updated
//-------------------------------------------------------------------------------------------------------------------------
   public void step(HistogramFrame xVelocityHistogram, HistogramFrame yVelocityHistogram, HistogramFrame speedHistogram) {
      verlet_solver.step();
      double totalKineticEnergy = 0;
      double v2 = 0;

      for(int i = 0; i<N; i++) {
	v2 = (state[(4 * i) + 1] * state[(4 * i) + 1]) + (state[(4 * i) + 3] * state[(4 * i) + 3]);

	if((v2 / 2) > (4 * initialKineticEnergy) ) {
		// if the velocity of the i-th particle is too high, and the user input bootsrap is true, then the particle
		// will be artificially slowed down to keep the system from blowing up
		highKE = true;
		if(bootstrap && ( (((v2/2)>(10 * lastKineticEnergy)) && (lastKineticEnergy != 0)) || ((v2/2) > (1000 * initialKineticEnergy))) ) {
			//System.out.println("Here");
			phi = Math.atan((state[4 * i + 1] / state[4 * i + 3]));
			if(state[4 * i + 1] < 0 && state[4 * i + 3] > 0) {
				phi = Math.PI - phi;
			} else if(state[4 * i + 1] < 0 && state[4 * i + 3] < 0) {
				phi = Math.PI + phi;
			} else if(state[4 * i + 1] > 0 && state[4 * i + 3] < 0) {
				phi = 2 * Math.PI - phi;
			}
			state[4 * i + 1] = Math.sqrt(initialKineticEnergy / N) * Math.cos(phi);
			state[4 * i + 3] = Math.sqrt(initialKineticEnergy / N) * Math.sin(phi); //bootstraps the kinetic energy down after an explosion
			resetAverages();
			steps = 1;		// get rid of the problems caused by an exploding temperature
		}
	}

	// calculate the totalk kinetic energy as well as append the current velocities and speed to the histograms
         totalKineticEnergy += (state[4*i+1]*state[4*i+1]+state[4*i+3]*state[4*i+3]);
         xVelocityHistogram.append(state[4*i+1]);
	 yVelocityHistogram.append(state[4 * i + 3]);
	 speedHistogram.append((state[4 * i + 1] * state[4 * i + 1] + state[4 * i + 3] * state[4 * i + 3]));
         state[4*i] = pbcPosition(state[4*i], Lx);
         state[4*i+2] = pbcPosition(state[4*i+2], Ly);
      }

      totalKineticEnergy *= 0.5; // since KE = m v^2 / 2 this finds the total kinetic energy on this step
      lastKineticEnergy = totalKineticEnergy;
      steps++;
      totalKineticEnergyAccumulator += totalKineticEnergy;
      totalKineticEnergySquaredAccumulator += totalKineticEnergy*totalKineticEnergy;
      t += dt;
   }

//draw() outputs the graphs for the LJParticles display
//--------------------------------------------------------------------------------------------------------------------------
   public void draw(DrawingPanel myWorld, Graphics g) {

      if(state==null) {
         return;
      }

      int pxRadius = Math.abs(myWorld.xToPix(radius)-myWorld.xToPix(0));
      int pyRadius = Math.abs(myWorld.yToPix(radius)-myWorld.yToPix(0));
      g.setColor(Color.red);
      for(int i = 0; i<N; i++) {	// this draws each particle
         int xpix = myWorld.xToPix(state[4*i])-pxRadius;
         int ypix = myWorld.yToPix(state[4*i+2])-pyRadius;
         g.fillOval(xpix, ypix, 2*pxRadius, 2*pyRadius);
      }


      if(isSilo) {			// this draws the aperture for the base of the silo
	int botLeftHole = myWorld.xToPix( ((Lx - holeDiam)/2) );
	int topHole = myWorld.yToPix(0.1);
	int xRange = myWorld.xToPix((Lx+holeDiam)/2) - myWorld.xToPix((Lx - holeDiam)/2);
	int yRange = myWorld.yToPix(0) - myWorld.yToPix(0.1);
	g.setColor(Color.BLUE);
	g.fillRect(botLeftHole, topHole, xRange , yRange);
      }

      g.setColor(Color.black);		// this draws the boundary for the container
      int xpix = myWorld.xToPix(0);
      int ypix = myWorld.yToPix(Ly);
      int lx = myWorld.xToPix(Lx)-myWorld.xToPix(0);
      int ly = myWorld.yToPix(0)-myWorld.yToPix(Ly);
      g.drawRect(xpix, ypix, lx, ly);
   }

//findNeighbors() notes all particles within 3 cut-off radii of each other
//-------------------------------------------------------------------------------------------------------------------------
  public void findNeighbors() {
		int i = 0;
	int j = 0;
	double dx = 0;
	double dy = 0;
	double rSq = 0;
	while(i < (N - 1) ) { // want to go only up to the second to last particle (i = (N-2)) so that
				// j being set to i + 1 doesn't step out of the array
		numberNeighbors[i] = 0;
		j = i + 1;
			while(j < N) {//this will stop at j = N - 1 the last particle in the array
				dx = pbcSeparation(state[4 * i] - state[4 * j], Lx);
				dy = pbcSeparation(state[(4 * i) + 2] - state[(4 * j) + 2], Ly);
				rSq = (dx * dx) + (dy * dy);
				if(rSq < (9 * rCutOffSq)) {//i.e. r < 3 rCutOff
					neighbors[i][numberNeighbors[i]] = j;// then j is i's neighbor
					numberNeighbors[i] = numberNeighbors[i] + 1;
				}
				j = j + 1;
			} //ends while(j < N)
		i = i + 1;
	} //ends while(i < (N - 1))
  } // ends findNeighbors()

}