1 --
    2 --      Patricia Fasel
    3 --      Los Alamos National Laboratory
    4 --      1990 August
    5 --
    6 module PushParticle (pushParticle) where
    7 
    8 import  PicType
    9 import  Consts
   10 import Array--1.3
   11 
   12 -- Phase IV: Particle push
   13 -- Each particle has an initial velocity determined by its motion during the
   14 --     previous timestep and is accelerated by the field induced by the 
   15 --     collection of particles
   16 --
   17 -- Compute the acceleration of each particle
   18 -- Find the maximum acceleration in x and y directions
   19 -- Determine a safe delta t, such that no particle will move a
   20 --     distance greater than the cell
   21 -- Compute new position and velocity of each particle
   22 
   23 pushParticle :: ParticleHeap -> Electric -> Value -> Value -> Value -> (Value, Value, ParticleHeap) 
   24 pushParticle ([], []) xyElec dt maxAcc maxVel = (maxAcc, maxVel, ([], []))
   25 pushParticle (((xPos,yPos):xyPos), ((xVel,yVel):xyVel))
   26              (xElec, yElec) dt maxAcc maxVel =
   27         (maxAcc'', maxVel'', 
   28             (((xPos',yPos'):xyPos'), ((xVel',yVel'):xyVel')))
   29         where
   30             i = truncate xPos
   31             j = truncate yPos
   32             i1 = (i+1) `rem` nCell
   33             j1 = (j+1) `rem` nCell
   34             dx = xPos - fromIntegral i
   35             dy = yPos - fromIntegral j
   36             xAcc = (charge/mass) * (xElec!(i,j)*(1-dy) + xElec!(i,j1)*dy)
   37             yAcc = (charge/mass) * (yElec!(i,j)*(1-dx) + yElec!(i1,j)*dx)
   38             xTV = xAcc*dt + xVel
   39             yTV = yAcc*dt + yVel
   40             xT = xTV*dt + xPos
   41             yT = yTV*dt + yPos
   42             maxAcc' = max maxAcc (max (abs xAcc) (abs yAcc))
   43             maxVel' = max maxVel (max (abs xTV) (abs yTV))
   44             (xVel',yVel') = (xTV, yTV)
   45             xPos' = if (xT >= fromIntegral nCell) 
   46                         then xT - fromIntegral nCell
   47                         else if (xT < 0.0) 
   48                             then xT + fromIntegral nCell
   49                             else xT
   50             yPos' = if (yT >= fromIntegral nCell) 
   51                         then yT - fromIntegral nCell
   52                         else if (yT < 0.0) 
   53                             then yT + fromIntegral nCell
   54                             else yT
   55             (maxAcc'', maxVel'', (xyPos', xyVel')) =
   56                 pushParticle (xyPos, xyVel) (xElec, yElec) dt maxAcc' maxVel'