package org.openpixi.pixi.physics.grid;

import org.openpixi.pixi.physics.Debug;
import org.openpixi.pixi.physics.particles.Particle;

/**
 *Interpolates current from the particles to the grid in a way s.t. the continuity equation 
 * divJ = - drho/dt
 *is fulfilled. This condition depends also on the way how the charge density is calculated.
 *This algorithm assumes area weighting used in the CloudInCell algorithm.
 *NOTE: On a coarser grid this algorithm will give a lower current compared to a finer grid
 *when the particle travels the same absolute distance in both cases. 
 */
public class ChargeConservingCIC extends CloudInCell {

	@Override
	public void interpolateToGrid(Particle p, Grid g, double tstep) {
		/**X index of local origin i.e. nearest grid point BEFORE particle push*/
		int xStart;
		/**Y index of local origin i.e. nearest grid point BEFORE particle push*/
		int yStart;
		/**X index of local origin i.e. nearest grid point AFTER particle push*/
		int xEnd;
		/**Y index of local origin i.e. nearest grid point AFTER particle push*/
		int yEnd;
		/**Normalized local x coordinate BEFORE particle push*/
		double x;
		/**Normalized local y coordinate BEFORE particle push*/
		double y;
		/**Normalized distance covered in X direction*/
		double deltaX;
		/**Normalized distance covered in X direction*/
		double deltaY;
		
        //we start by shifting the particle by (-0.5, -0.5)
        //since the grid for the current is shifted by (0.5,0.5)
        p.addPrevX(-0.5*g.getCellWidth());
        p.addPrevY(-0.5*g.getCellHeight());
        p.addX(-0.5*g.getCellWidth());
        p.addY(-0.5*g.getCellHeight());
		
		x = p.getPrevX() / g.getCellWidth();
		y = p.getPrevY() / g.getCellHeight();
		
		xStart = (int) Math.floor(x + 0.5);
		yStart = (int) Math.floor(y + 0.5);
		
		deltaX = p.getX() / g.getCellWidth();
		deltaY = p.getY() / g.getCellHeight();
		
		xEnd = (int) Math.floor(deltaX + 0.5);
		yEnd = (int) Math.floor(deltaY + 0.5);
		
		deltaX -= x;
		deltaY -= y;
		
		x -= xStart;
		y -= yStart;
		
		//check if particle moves further than one cell
		if (Debug.asserts) {
			assert (Math.abs(deltaX) <= g.getCellWidth()) & (Math.abs(deltaY) <= g.getCellHeight()): "particle too fast";
		}

		//4-boundary move?
		if (xStart == xEnd && yStart == yEnd) {
			fourBoundaryMove(xStart, yStart, x, y, deltaX, deltaY, p, g, tstep);
			}
		//7-boundary move?
		else if (xStart == xEnd || yStart == yEnd) {
				sevenBoundaryMove(x, y, xStart, yStart, xEnd, yEnd, deltaX, deltaY, p, g, tstep);
			}
			// 10-boundary move
				else {
					tenBoundaryMove(x, y, xStart, yStart, xEnd, yEnd, deltaX, deltaY, p, g, tstep);
				}
        p.addPrevX(0.5*g.getCellWidth());
        p.addX(0.5*g.getCellWidth());
        p.addPrevY(0.5*g.getCellHeight());
        p.addY(0.5*g.getCellHeight());
	}

	/**
	 * Writes currents Jx and Jy to grid.
	 * @param lx x index of local origin
	 * @param ly y index of local origin
	 * @param x local x coordinate relative to lx BEFORE particle push
	 * @param y local y coordinate relative to ly BEFORE particle push
	 * @param deltaX x distance covered by particle (not absolute but only for this 4-boundary move)
	 * @param deltaY y distance covered by particle (not absolute but only for this 4-boundary move)
	 * @param tstep
	 */
	private void fourBoundaryMove(int lx, int ly, double x, double y,
								  double deltaX, double deltaY, Particle p, Grid g, double tstep) {
		//A few cancellations were made to reduce computation time. Till this point the algorithm has
		//calculated the area that swept over a cell boundary for a normalized grid (i.e. unit square cells).
		//and unit square charges. This area needs to be denormalized and then multiplied with the charge
		//density. But these operations cancel and no further calculations need to be done.
		//g.addJx(lx, 	ly - 1, p.getCharge() * deltaX * ((1 - deltaY) / 2 - y));
		//g.addJx(lx, 	ly, 	p.getCharge() * deltaX * ((1 + deltaY) / 2 + y));
		//g.addJy(lx - 1, ly, 	p.getCharge() * deltaY * ((1 - deltaX) / 2 - x));
		//g.addJy(lx, 	ly, 	p.getCharge() * deltaY * ((1 + deltaX) / 2 + x));
		/*
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), p.getCharge() * deltaX * ((1 - deltaY) / 2 - y) * g.getCellWidth() / tstep);
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + 1 + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge()* deltaX * ((1 + deltaY) / 2 + y) * g.getCellWidth() / tstep);
		g.addJy((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 - deltaX) / 2 - x) * g.getCellHeight() / tstep);
		g.addJy((lx +1 + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 + deltaX) / 2 + x) * g.getCellHeight() / tstep);
		*/
		
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), p.getCharge() * deltaX * ((1 - deltaY) / 2 - y) / g.getCellWidth() / tstep);
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + 1 + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge()* deltaX * ((1 + deltaY) / 2 + y) / g.getCellWidth() / tstep);
		g.addJy((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 - deltaX) / 2 - x) / g.getCellHeight() / tstep);
		g.addJy((lx +1 + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 + deltaX) / 2 + x) / g.getCellHeight() / tstep);
		
		/*
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), p.getCharge() * deltaX * ((1 - deltaY) / 2 - y));
		g.addJx((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + 1 + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge()* deltaX * ((1 + deltaY) / 2 + y));
		g.addJy((lx + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 - deltaX) / 2 - x));
		g.addJy((lx +1 + g.getNumCellsX())%g.getNumCellsX(), (ly + g.getNumCellsY())%g.getNumCellsY(), 	p.getCharge() * deltaY * ((1 + deltaX) / 2 + x));
		*/

	}

	private void sevenBoundaryMove(double x, double y, int xStart, int yStart, int xEnd, int yEnd,
								   double deltaX, double deltaY, Particle p, Grid g, double tstep) {
		//7-boundary move with equal y?
		if (yStart == yEnd) {
			//particle moves right?
			if (xEnd > xStart) {

				double deltaX1 = 0.5 - x;
				double deltaY1 = (deltaY / deltaX) * deltaX1;
				fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

				deltaX -= deltaX1;
				deltaY -= deltaY1;
				y += deltaY1;
				fourBoundaryMove(xEnd, yEnd, - 0.5, y, deltaX, deltaY, p, g, tstep);
				
			}
			//particle moves left
			else {

				double deltaX1 = -(0.5 + x);
				double deltaY1 = (deltaY / deltaX) * deltaX1;
				fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

				deltaX -= deltaX1;
				deltaY -= deltaY1;
				y += deltaY1;
				fourBoundaryMove(xEnd, yEnd, 0.5, y, deltaX, deltaY, p, g, tstep);

			}
		}
		//7-boundary move with equal x?
		if (xStart == xEnd) {
			//particle moves up?
			if (yEnd > yStart) {

				double deltaY1 = 0.5 - y;
				double deltaX1 = deltaX  * (deltaY1 / deltaY);
				fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

				deltaX -= deltaX1;
				deltaY -= deltaY1;
				y += deltaY1;
				fourBoundaryMove(xEnd, yEnd, x, -0.5, deltaX, deltaY, p, g, tstep);

			}
			//particle moves down
			else {

				double deltaY1 = -(0.5 + y);
				double deltaX1 = (deltaX / deltaY) * deltaY1;
				fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

				deltaX -= deltaX1;
				deltaY -= deltaY1;
				y += deltaY1;
				fourBoundaryMove(xEnd, yEnd, x, 0.5, deltaX, deltaY, p, g, tstep);

			}
		}

	}

	private void tenBoundaryMove(double x, double y, int xStart, int yStart, int xEnd, int yEnd,
								 double deltaX, double deltaY, Particle p, Grid g, double tstep) {
		//moved right?
		if (xEnd == (xStart+1)) {
			//moved up?
			if (yEnd == (yStart+1)) {

				double deltaX1 = 0.5 - x;

				//lower local origin
				if(((deltaY / deltaX) * deltaX1 + y) < 0.5) {

					double deltaY1 = (deltaY / deltaX) * deltaX1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaY2 = 0.5 - y - deltaY1;
					double deltaX2 = (deltaX1 / deltaY1) * deltaY2;
					y += deltaY1;
					fourBoundaryMove(xStart+1, yStart, -0.5, y, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					x = deltaX2 - 0.5;
					fourBoundaryMove(xEnd, yEnd, x, -0.5, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 >= 0: deltaX1;
						assert deltaY1 >= 0: deltaY1;
						assert y >= 0 && y <= 0.5;
						assert deltaX2 >= 0: deltaX2;
						assert deltaY2 >= 0: deltaY2;
						assert deltaX >= 0: deltaX;
						assert deltaY >=0: deltaY;
						assert x <= 0 && x >= -0.5;
					}
				}
				//upper local origin
				else {

					double deltaY1 = 0.5 - y;
					deltaX1 = (deltaX / deltaY) * deltaY1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaX2 = 0.5 - x - deltaX1;
					double deltaY2 = (deltaY1 / deltaX1) * deltaX2;
					x += deltaX1;
					fourBoundaryMove(xStart, yStart+1, x, -0.5, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					y = deltaY2 - 0.5;
					fourBoundaryMove(xEnd, yEnd, -0.5, y, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 >= 0: deltaX1;
						assert deltaY1 >= 0: deltaY1;
						assert x >= 0 && x <= 0.5;
						assert deltaX2 >= 0: deltaX2;
						assert deltaY2 >= 0: deltaY2;
						assert deltaX >=0: deltaX;
						assert deltaY >=0: deltaY;
						assert y <= 0 && y >= -0.5;
					}
				
				}
			}
			//moved down
			else {

				double deltaY1 = -(0.5 + y);

				//lower local origin
				if(((deltaX / deltaY) * deltaY1 + x) < 0.5) {

					double deltaX1 = (deltaX / deltaY) * deltaY1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaX2 = 0.5 - x - deltaX1;
					double deltaY2 = (deltaY / deltaX) * deltaX2;
					x += deltaX1;
					fourBoundaryMove(xStart, yStart-1, x, 0.5, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					y = 0.5 + deltaY2;
					fourBoundaryMove(xEnd, yEnd, -0.5, y, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaY1 <= 0: deltaY1;
						assert x >= 0 && x <= 0.5;
						assert deltaX1 >= 0: deltaX1;
						assert deltaX2 >= 0: deltaX2;
						assert deltaY2 <= 0: deltaY2;
						assert deltaX >= 0: deltaX;
						assert deltaY <=0: deltaY;
						assert y >= 0 && y <= 0.5;
					}

				}
				//upper local origin
				else {

					double deltaX1 = 0.5 - x;
					deltaY1 = (deltaY / deltaX) * deltaX1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaY2 = -(0.5 + y + deltaY1);
					double deltaX2 = (deltaX1 / deltaY1) * deltaY2;
					y += deltaY1;
					fourBoundaryMove(xStart+1, yStart, -0.5, y, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					x = deltaX2 - 0.5;
					fourBoundaryMove(xEnd, yEnd, x, 0.5, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 >= 0: deltaX1;
						assert deltaY1 <= 0: deltaY1;
						assert y <= 0 && y >= -0.5;
						assert deltaX2 >= 0: deltaX2;
						assert deltaY2 <= 0: deltaY2;
						assert deltaX >= 0: deltaX;
						assert deltaY <=0: deltaY;
						assert x <= 0 && x >= -0.5;
					}
				}
			}

		}
		//moved left
		else {
			//moved up?
			if (yEnd == (yStart+1)) {

				double deltaX1 = -(0.5 + x);
				//lower local origin
				if(((deltaY / deltaX) * deltaX1 + y) < 0.5) {

					double deltaY1 = (deltaY / deltaX) * deltaX1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaY2 = 0.5 - y - deltaY1;
					double deltaX2 = (deltaX1 / deltaY1) * deltaY2;
					y += deltaY1;
					fourBoundaryMove(xStart-1, yStart, 0.5, y, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					x = 0.5 + deltaX2;
					fourBoundaryMove(xEnd, yEnd, x, -0.5, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 <= 0: deltaX1;
						assert deltaY1 >= 0: deltaY1;
						assert y >= 0 && y <= 0.5;
						assert deltaX2 <= 0: deltaX2;
						assert deltaY2 >= 0: deltaY2;
						assert deltaX <= 0: deltaX;
						assert deltaY >=0: deltaY;
						assert x >= 0 && x <= 0.5;
					}
				}
				//upper local origin
				else {

					double deltaY1 = 0.5 - y;
					deltaX1 = (deltaX / deltaY) * deltaY1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaX2 = -(0.5 + x + deltaX1);
					double deltaY2 = (deltaY1 / deltaX1) * deltaX2;
					x += deltaX1;
					fourBoundaryMove(xStart, yStart+1, x, -0.5, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					y = deltaY2 - 0.5;
					fourBoundaryMove(xEnd, yEnd, 0.5, y, deltaX, deltaY,p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 <= 0: deltaX1;
						assert deltaY1 >= 0: deltaY1;
						assert x <= 0 && x >= -0.5;
						assert deltaX2 <= 0: deltaX2;
						assert deltaY2 >= 0: deltaY2;
						assert deltaX <=0: deltaX;
						assert deltaY >=0: deltaY;
						assert y <= 0 && y >= -0.5;
					}
				}
			}
			//moved down
			else {

				double deltaY1 = -(0.5 + y);
				//lower local origin
				if((-(deltaX / deltaY) * deltaY1 - x) < 0.5) {

					double deltaX1 = (deltaX / deltaY) * deltaY1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1,p, g, tstep);

					double deltaX2 = -(0.5 + x + deltaX1);
					double deltaY2 = (deltaY / deltaX) * deltaX2;
					x += deltaX1;
					fourBoundaryMove(xStart, yStart-1, x, 0.5, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					y = 0.5 + deltaY2;
					fourBoundaryMove(xEnd, yEnd, 0.5, y, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaY1 <= 0: deltaY1;
						assert deltaX1 <= 0: deltaX1;
						assert x <= 0 && x >= -0.5;
						assert deltaX2 <= 0: deltaX2;
						assert deltaY2 <= 0: deltaY2;
						assert deltaX <= 0: deltaX;
						assert deltaY <= 0: deltaY;
						assert y >= 0 && y <= 0.5;
					}

				}
				//upper local origin
				else {

					double deltaX1 = -(0.5 + x);
					deltaY1 = (deltaY / deltaX) * deltaX1;
					fourBoundaryMove(xStart, yStart, x, y, deltaX1, deltaY1, p, g, tstep);

					double deltaY2 = -(0.5 + y + deltaY1);
					double deltaX2 = (deltaX1 / deltaY1) * deltaY2;//Error source!!!
					y += deltaY1;
					fourBoundaryMove(xStart+1, yStart, 0.5, y, deltaX2, deltaY2, p, g, tstep);

					deltaX -= (deltaX1 + deltaX2);
					deltaY -= (deltaY1 + deltaY2);
					x = 0.5 + deltaX2;
					fourBoundaryMove(xEnd, yEnd, x, 0.5, deltaX, deltaY, p, g, tstep);

					if (Debug.asserts) {
						assert deltaX1 <= 0: deltaX1;
						assert deltaY1 <= 0: deltaY1;
						assert y <= 0 && y >= -0.5;
						assert deltaX2 <= 0: deltaX2;
						assert deltaY2 <= 0: deltaY2;
						assert deltaX <= 0: deltaX;
						assert deltaY <= 0: deltaY;
						assert x >= 0 && x <= 0.5;
					}
				}
			}
		}
	}
}
