-- 
--      Patricia Fasel
--      Los Alamos National Laboratory
--      1990 August
--
module ChargeDensity (chargeDensity) where

import	PicType
import	Consts
import Array--1.3

-- Phase I: calculate the charge density rho
-- Each particle represents an amount of charge distributed over an entire cell.
-- In discretizing the charge density at the grid points we split the cell into
--    four rectangles and assign to each corner an amount of charge proportional
--    to the area of the opposite diagonal sub rectangle
-- So each particle contributes a portion of charge to four quadrants
-- particle position is (i+dx, j+dy)
-- nw quadrant (1-dx)(1-dy)(charge) is added to rho(i,j)
-- ne quadrant (dx)(1-dy)(charge) is added to rho(i,j+1)
-- sw quadrant (1-dx)(dy)(charge) is added to rho(i+1,j)
-- se quadrant (dx)(dy)(charge) is added to rho(i+1,j+1)
-- wrap around used for edges and corners

chargeDensity :: ParticleHeap -> Rho
chargeDensity (xyPos, xyVel) =
	accumArray (+) 0 ((0,0), (n,n)) (accumCharge xyPos)
	where
	    n = nCell-1


-- for each particle, calculate the distribution of density
-- based on (x,y), a proportion of the charge goes to each rho

accumCharge :: [Position] -> [MeshAssoc]
accumCharge [] = []
accumCharge ((x,y):xys) =
	[((i ,j ) , charge * (1-dx) * (1-dy))] ++
	[((i',j ) , charge * dx * (1-dy))] ++
	[((i ,j') , charge * (1-dx) * dy)] ++
	[((i',j') , charge * dx * dy)] ++
	accumCharge xys
	where
	    i = truncate x
	    i' = (i+1) `rem` nCell
	    j = truncate y
	    j' = (j+1) `rem` nCell
	    dx = x - fromIntegral i
	    dy = y - fromIntegral j