ASTR 5540 Math Meth Fall 2008. Problem Set 9. Due Mon Nov 3
1. 1-D PM code
In the language of your choice (IDL, Mathematica, c, ...), write a 1-dimensional gravita-
tional Particle-Mesh (PM) code. Choose whatever units of length, mass, and time make
programming easiest.
(a) Functions
Your code should contain functions that implement the following:
1. A cloud-in-cell (CIC) smoothing window W(x);
2. Wavenumber of a Fourier mode;
3. Mesh density given particle positions;
4. FFT of mesh density;
5. Potential given density, in Fourier space;
6. Acceleration given potential, in Fourier space;
7. FFT acceleration back on to mesh;
8. Acceleration at particle positions given mesh acceleration.
[Hints: Function 2 depends on how the FFT you use stores the FT. In mathematica for exam-
ple, the FT is stored as a complex array {˜a0,˜a1,˜a2, ..., ˜a[N
2],˜a[N
2]+1−N, ..., ˜a−2,˜a−1}, so that the
wavenumbers are {0,1,2, ..., [N
2],[N
2]+1−N, ..., −2,−1}. FFTs that input a real array typi-
cally return {˜a0,Re ˜a1,Im ˜a1, ..., Re ˜a[N
2]}, in which case the wavenumbers are {0,1,1, ..., [N
2]}.
Function 5 should set the zero’th Fourier mode of the potential to zero, ˜
φ0= 0. Function
6 should set the zero’th and [N
2]’th (for even N) Fourier mode of the acceleration to zero,
˜g0= ˜g[N
2]= 0. The sign of the acceleration ˜gin Fourier space depends on the phase con-
vention of the FFT; for example mathematica’s phase convention is minus the convention
adopted in class.]
(b) Leap-frog integrator
Implement a leap-frog integrator for the positions xand velocities vof particles:
vi+1
2=vi−
1
2+gi∆t ,
xi+1 =xi+vi+1
2∆t , (1.1)
where ∆tis a suitable time step.
(c) Integrate and plot
Choose a reasonable mesh-size N. Choose a small number of particles, say 2 to 4. Give
the particles random initial positions, and zero initial velocities. Choose a time step ∆t.
Integrate for a sufficient number of steps, at least 100. Plot the positions of the particles as
a function of time t.
1