We describe a method for simulating first-order reversal curve (FORC) diagrams of interacting single-domain particles. Magnetostatic interactions are calculated in real space, allowing simulations to be performed for particle ensembles with arbitrary geometry. For weakly interacting uniaxial particles, the equilibrium magnetization at each field step is obtained by direct solution of the Stoner-Wohlfarth model, assuming a quasi-static distribution of interaction fields. For all other cases, the equilibrium magnetization is calculated using an approximate iterated solution to the Landau-Lifshitz-Gilbert equation. Multithreading is employed to allow multiple curves to be computed simultaneously, enabling FORC diagrams to be simulated in reasonable time using a standard desktop computer. Statistical averaging and post processing lead to simulated FORC diagrams that are comparable to their experimental counterparts. The method is applied to several geometries of relevance to rock and environmental magnetism, including densely packed random clusters and partially collapsed chains. The method forms the basis of FORCulator, a freely available software tool with graphical user interface that will enable FORC simulations to become a routine part of rock magnetic studies.