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

/**
 * SiloParticles extends LJParticles to model granular particles within a silo with a hole in the base
 *
 * Daniel Costantino
 */
public class SiloParticles extends LJParticles {

	public SiloParticles() {
		super();
	}

//computeAcceleration is slightly modified since the boundary conditions are different from LJParticles, but most of it
//is the same as in LJParticles-------------------------------------------------------------------------------------------- 
   public void computeAcceleration() {
      int j = 0;
      int k = 0;
      for(int i = 0; i<N; i++) {
         ax[i] = 0;
         ay[i] = 0;
      } // initializes acceleration
      int randParticle = (int)(N * Math.random());

      for(int i = 0; i<=N-1; i++) {
	double dx = 0;
	double dy = 0;
	double r2 = 0;
	double r = 0;
	double fx = 0;
	double fy = 0;
	double friction = 0;
	double vx = 0;
	double vy = 0;
	double oneOverR2 = 0;
	double oneOverR6 = 0;
	double oneOverRCO2 = 1.0 / rCutOffSq;
	double oneOverRCO6 = oneOverRCO2 * oneOverRCO2 * oneOverRCO2;
	double v2 = 0;
	double dampingV2 = 0;
	double fOverR = 0;
	double fatRC = 0;
	double bootStrapV = 0;

	//smallR = false;

         while(k < numberNeighbors[i] && (consider[i] == 1) )  {

	    j = neighbors[i][k];

            dx = state[4*i]-state[4*j];		//in a silo, the separation between grains is just their separation on
            dy = state[4*i+2]-state[4*j+2];	//the grid, no periodic boundary conditions apply

            r2 = dx*dx+dy*dy;
	    r = Math.sqrt(r2);
	    oneOverR2 = 1.0/r2;
	    oneOverR6 = oneOverR2*oneOverR2*oneOverR2;

	    if( r2 > rCutOffSq || (nuclear && ((dx < 0.15) && (dy < 0.15)) ) ) {
			fx = 0;
			fy = 0;
			friction = 0;
	    } else if((r2 < rCutOffSq) && (consider[j] == 1) ) {
		fOverR = 48.0*oneOverR6*(oneOverR6-0.5)*oneOverR2;
		fatRC = 48 * oneOverRCO6 * (oneOverRCO6 - 1/2) * oneOverRCO2;
		fx = (fOverR - fatRC) * dx;
		fy = (fOverR - fatRC) * dy;
		vx = state[(4 * i) + 1] - state[(4 * j) + 1];
		vy = state[(4 * i) + 3] - state[(4 * j) + 3];
		friction = -gamma * (vx * dx + vy * dy) / r2;
		    totalPotentialEnergyAccumulator += 4.0*((oneOverR6*oneOverR6-oneOverR6) - (oneOverRCO6 * oneOverRCO6 - oneOverRCO6));
							 // to avoid discontinuity in energy

	    }


	    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));
	    virialAccumulator += dx*fx+dy*fy;
	    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) {
		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) * Math.cos(theta);
		state[4 * i + 3] = Math.sqrt(2.0 * desiredKE) * Math.sin(theta);
	}

	ay[i] = ay[i] - g;	// the particles feel a downward force due to gravity

      } // ends for(i < N - 1)
   }

//siloPositionX takes care of the boundary conditions on either side of the container
//------------------------------------------------------------------------------------------------------------------------
   private double[] siloPositionX(double xPos, double xVel, double xLength) {
	if(xPos < 0) {
		xPos = -xPos;
		while(xPos > xLength) {
			xPos = xPos - xLength;
		}
		if( (xPos < (xLength / 2)) && (xVel < 0) ) {
				//the particles should go through an elastic collision with the walls
			xVel = -xVel;
		}
	} else if(xPos > xLength) {
		xPos = 2 * xLength - xPos;
		while(xPos < 0) {
			xPos = xLength + xPos;
		}
		if( (xPos > (xLength / 2)) && (xVel > 0) ) {
			xVel = -xVel;
		}
	}
	double xInfo[] = new double[2];	//new xPosition and xVelocity are passed back to step()
	xInfo[0] = xPos;
	xInfo[1] = xVel;
	return xInfo;
   }

//siloPositionY takes care of the boundary conditions at the top and bottom of the container and removes particles if they
//pass through the hole----------------------------------------------------------------------------------------------------
   private double[] siloPositionY(double xPos, double yPos, double xVel, double yVel, double holeSize, double xLength, double yLength, int numPart) {
	if(yPos < 0) {
		if( (xPos < (xLength - holeSize) / 2) || (xPos > (xLength + holeSize) / 2) ) {
			yPos = -yPos;
			while(yPos > yLength) {
				yPos = yPos - yLength;
			}
			if( (yPos < (yLength / 2)) && (yVel < 0)) {
				yVel = -yVel;
			}
		} else {
				// if a particle passes through the hole, it is dropped from the ceiling directly above where
				// it passed through the hole and starts from rest
			yPos = yLength + yPos;
			while(yPos < 0) {
				yPos = yPos - yLength;
			}
			
			yVel = 0;
			xVel = 0;

			numPart = numPart - 1;
		}
	} else if(yPos > yLength) {
		yPos = 2 * yLength - yPos;
		while(yPos < 0) {
			yPos = yLength + yPos;
		}
		if((yPos > (yLength / 2)) && yVel > 0) {
			yVel = -yVel;
		}
	}
	double yInfo[] = new double[4];	//new yPosition, xVelocity, yVelocity, and number of particles are passed back to
	yInfo[0] = yPos;		//step()
	yInfo[1] = yVel;
	yInfo[2] = numPart;
	yInfo[3] = xVel;
	return yInfo;
   }


//step goes through a time step, calls the Verlet solver and deals with boundary conditions-------------------------------
//------------------------------------------------------------------------------------------------------------------------
   public void step(HistogramFrame xVelocityHistogram, HistogramFrame yVelocityHistogram, HistogramFrame speedHistogram) {
      verlet_solver.step();
      double totalKineticEnergy = 0;
      double v2 = 0;
      double xInformation[] = new double[2];
      double yInformation[] = new double[3];

      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
		}
	 }
 
         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]));

	 //call up the boundary conditions and get the velocities and positions dealt with correctly: 
         xInformation = siloPositionX(state[4*i], state[4 * i + 1], Lx);
	 state[4 * i] = xInformation[0];
	 state[4 * i + 1] = xInformation[1];


	 yInformation = siloPositionY(state[4 * i], state[4*i+2], state[4 * i + 1], state[4 * i + 3], holeDiam, Lx, Ly, nPart);

	 state[4 * i + 2] = yInformation[0];
	 state[4 * i + 3] = yInformation[1];
	 nPart = (int) yInformation[2];
	 state[4 * i + 1] = yInformation[3];

      }

      totalKineticEnergy *= 0.5;
      lastKineticEnergy = totalKineticEnergy;
      steps++;
      totalKineticEnergyAccumulator += totalKineticEnergy;
      totalKineticEnergySquaredAccumulator += totalKineticEnergy*totalKineticEnergy;
      t += dt;
   }
}//ends SiloParticles