there is no particular need to have separate phantomsetup utilities for 1, 2 and many stars
It would be better to combine setup_star.f90, setup_binary.f90 and setup_hierarchical.f90 into a single all-singing-all-dancing setup for 1, 2 or many stars, with any combination of sink particles and stars setup from stellar profiles.
In principle this is not very much work, and would make a single well-supported setup for stellar interaction simulations. The work required is:
shift any remaining functionality from setup_star.f90 into setup_binary.f90 (mainly read of equation of state options from .setup file)
add an option in the .setup file specifying the number of stars to setup (1, 2 or >=3)
handle the n=1 case in setup_binary.f90
extend setup_binary to call set_hierarchical instead of set_binary to handle the >= 3 case
rename setup_binary to setup_stars.f90 as the mother-of-all star setup routines
it may also be possible to delete SETUP=common and SETUP=polytrope (i.e. make them aliases to SETUP=star). These don't seem to add any functionality that is not already present in the other routines
there is no particular need to have separate phantomsetup utilities for 1, 2 and many stars
It would be better to combine setup_star.f90, setup_binary.f90 and setup_hierarchical.f90 into a single all-singing-all-dancing setup for 1, 2 or many stars, with any combination of sink particles and stars setup from stellar profiles.
In principle this is not very much work, and would make a single well-supported setup for stellar interaction simulations. The work required is: