This set of MATLAB routines simulates a three-dimensional hard sphere packing with periodic boundary conditions. It initializes a uniform lattice, and then uses the Metropolis algorithm to anneal the particle locations for many iterations. It then calculates $g_2(r)$, or the pairwise radial distribution function, averaged over many centers in the set. For large initial particle densities, the resulting $g_2$ matches the analytic approximation given by the Percus-Yevick approximation to arbitrary precision given a high enough particle density.