Determine position of the proper "dynamical center of mass" of a system, e.g. using the 50% most bound particles. For this, it's best to save the potentials of the particles during the integration. Most N-body integrators will have an option to let you save the potentials of the particles. In hackcode1 one has to use the options=mass,phase,phi keyword to save this relevant information.
Here is an example of analyzing the resulting output dataset of the merger we ran earlier. Using snapcenter the snapshot is centered weighing each particle by the third power of the potential, followed by a computation of the lagrangian radii using snapmradii. The resulting table is piped into tabmath and the logarithm of all radii is computed (the first columns contains the time) after which the converted table is piped into tabplot.
set rad=0.01,0.05:0.95:0.05,0.99
set nrad=(`nemoinp $rad | wc -w`)
@ nrad++
snapcenter r001.dat - '-phi*phi*phi' report=f |\
snapmradii - $rad |\
tabmath - - %1,`nemoinp 2:$nrad format="log(%%%.f)" separ=,` all |\
tabplot - 1 2:$nrad line=1,1 yapp=2/xs