The current implementation of the reionization algorithm consists of one very long function (1500 lines). There is a lot of similar computation separated by flags (USE_MASS_DEPENDENT_ZETA, USE_MINI_HALOS, USE_HALO_FIELD) so there is a lot of opportunity to make this function more modular.

Functions could be created for initialising the required grids, finding maxima/minima for table generation, calculating ionizing emissivity from filtered grids at one radius, calculating recombinations etc.

By refactoring and removing redundant code here it will become much easier to add new models and detect issues in the future.