One dimensional Particle-Mesh (PM) code (template)

Define functions

(* define Cloud-in-Cell window *)
Clear[window];
window[x_]:=

(* define function giving wavenumber of Fourier mode *)
Clear[wavenumber];
wavenumber[m_]:=

(* define function giving density on mesh given a list of positions *)
Clear[meshdensity];
meshdensity[position_]:=

(* define function giving Fourier transform of mean-subtracted density; setting zero'th mode to zero is equivalent to subtracting mean density *)
Clear[meshdensityk];
meshdensityk[meshdensity_]:=

(* define function giving potential given mean-subtracted density in Fourier space; the zero'th mode of the potential is set to zero *)
Clear[meshpotentialk];
meshpotentialk[meshdensityk_]:=

(* define function giving acceleration given potential in Fourier space; (N/2)'th mode of acceleration is set to zero to ensure that acceleration is real *)
Clear[meshaccelerationk];
meshaccelerationk[meshpotentialk_]:=

(* define function that gives acceleration on mesh *)
Clear[meshacceleration];
meshacceleration[meshaccelerationk_]:=InverseFourier[meshaccelerationk]

(* define function that gives accelerations at positions of particles, given acceleration on mesh *)
Clear[positionacceleration];
positionacceleration[position_,meshacceleration_]:=

(* define function that gives accelerations at positions of particles, given positions of particles *)
Clear[acceleration];
acceleration[position_]:=positionacceleration[position,meshacceleration[meshaccelerationk[meshpotentialk[meshdensityk[meshdensity[position]]]]]]

Leap-frog integrator

Clear[leapfrog];
SetAttributes[leapfrog,HoldAll];
leapfrog[position_,velocity_,acceleration_,Δt_]:=(
position+=velocity Δt;
velocity+=acceleration Δt;
)

Set initial conditions and integrate

(* choose number of mesh points *)
Clear[Nmesh];
Nmesh=(2^8)

(* choose number of particles *)
Clear[Nparticle];
Nparticle=4

(* list of initial positions drawn uniformly in [0,Nmesh] *)
Clear[initialposition];
initialposition=Table[RandomReal[Nmesh],{i,1,Nparticle}]

(* list of initial velocities set to zero *)
Clear[initialvelocity];
initialvelocity=Table[0,{i,1,Nparticle}]

(* integrate positions and velocities *)
Clear[Δt,nstep];
Δt=4;
nstep=200;
Clear[position,positions,velocity,velocities];
position=initialposition;
velocity=initialvelocity;
positions=position;
velocities=velocity;
For[istep=1,istep≤nstep,istep++,
(* call leap-frog integrator *)
leapfrog[position,velocity,acceleration[position],Δt];
(* update lists of positions and velocities *)
Do[positions[[i]]=Flatten[{positions[[i]],position[[i]]}],{i,1,Nparticle}];
Do[velocities[[i]]=Flatten[{velocities[[i]],velocity[[i]]}],{i,1,Nparticle}];
]

ListLinePlot[positions,PlotRange->{0,Nmesh}]

ListLinePlot[velocities,PlotRange->All]