This code carries out the spherical-collapse calculation for standard cosmological models as well as for dark energy models when the dark energy can be taken to be spatially homogeneous. The calculation is valid on sub-horizon scales and takes a top-hat perturbation to exist in an otherwise featureless cosmos and follows its evolution into the non-linear regime where it reaches a maximum size and then recollapses. The code provides the user with the linear-collapse threshold (delta_c) and the virial overdensity (Delta_v) for the collapsed halo over a range of cosmic scale factors.
The user specifies a cosmological model and the code then solves both the linear and non-linear differential equations for the overdensity of the top-hat (using an RK4 method). The calculation stops when the overdensity of the top hat in the non-linear calculation becomes infinite and the value of delta_c is read off from the linear. The radius of the virialised halo is taken to be half of the radius at turn around, in accordance with the virial theorem in the Einstein-de Sitter case; the virial overdensity is then calculated at the collapse time. This is different from some other authors who include the cosmological constant in the virial condition, but it wasn't obvious to me how to extend this to dark energy models so I settled for the simple option of keeping the Einstein-de-Sitter virial radius condition. When including Lambda, for Omega_m=0.3, the virial density calculated is around 10 per cent higher than that provided by my code because the halo has to collapse further to counter the outward gravitational push of Lambda.
The code provides the user with dcDV.dat; a data file containing the scale factor at the time of collapse - a, matter density - Omega_m(a), delta_c(a), Delta_v(a), linear growth - g(a), logarithmic growth rate - f(a) and the integrated growth - G(a). These are provided for approximately 200 values of a. The user does not specify these a values and instead should interpolate from the data file to their desired a value. The integrations are accurate to 4 significant figures for delta_c (1.686 is recovered for EdS) and 3 significant figures for Delta_v (177 is recovered for EdS) and are slightly biased in the next significant figure, there is also some numerical noise that is particularly obvious in the delta_c calculation, but well below the level relevant for my purposes.
There are certainly cleverer ways to do this calculation; for example, using the top-hat radius rather than overdensity as the dynamical variable or using adaptive time-stepping in the ODE solver. However, my simplistic first attempt at the calculation was accurate enough for my purposes so I didn't invest any effort in making the calculation better. I am able to get the Einstein-de-Sitter calculation to run in ~20 seconds on my 2015 Macbook. I'm sure this runtime could be substantially improved if someone had some free afternoons.
The code compiles with gfortran, and probably with ifort too, and does not need to be linked to any libraries, so you can just gfortran collapse.f90 and off you go. Please let me know of any problems as soon as possible.