Script to compare Clustertools analysis to internal AMUSE analysis

compare_clustertools_to_amuse.py

"""
Script to compare Clustertools analysis to native AMUSE analysis
"""
# requires amuse-framework, amuse-masc and clustertool to be installed

import numpy as np
import matplotlib.pyplot as plt
import clustertools as ctools
from amuse.units import units, nbody_system
from amuse.units.quantities import VectorQuantity
from amuse.ext.masc import new_star_cluster

np.random.seed(17082022)


def main():
    mass_unit = units.MSun
    converter = nbody_system.nbody_to_si(1 | units.MSun, 1 | units.pc)
    stars = new_star_cluster()
    cluster = ctools.StarCluster(ctype='amuse')
    cluster.add_stars(
        stars.x, stars.y, stars.z,
        stars.vx, stars.vy, stars.vz,
        stars.mass,
        stars.key,
    )
    cluster.analyze(sortstars=True)
    results = {}
    results["Velocity"] = (
        VectorQuantity.new_from_scalar_quantities(*cluster.v)
        / stars.velocity.lengths()
    ).mean()
    results["Total mass"] = (cluster.mtot | mass_unit) / stars.total_mass()
    results["Mean mass"] = (cluster.mmean | mass_unit) / stars.mass.mean()
    results["Maximum radius"] = cluster.rmax / stars.total_radius()
    results["Mean radius"] = cluster.rmean / stars.virial_radius()
    stars.distance = stars.position.lengths()
    stars = stars.sorted_by_attribute("distance")
    cumulative_mass = stars.mass.cumsum()
    d10 = stars[
        cumulative_mass < 0.1 * stars.total_mass()
    ].distance.max()
    d50 = stars[
        cumulative_mass < 0.5 * stars.total_mass()
    ].distance.max()
    results["Half mass radius (vs AMUSE)"] = cluster.rm / stars.LagrangianRadii(
        unit_converter=converter, mf=[0.5,]
    )[0][0]
    results["Half mass radius (vs manual)"] = cluster.rm / d50
    results["10% Lagrange radius (vs AMUSE)"] = cluster.rm / stars.LagrangianRadii(
        unit_converter=converter, mf=[0.1,]
    )[0][0]
    results["10% Lagrange radius (vs manual)"] = cluster.rm / d10

    for key in results.keys():
        print(f"{key}: {results[key]}")

    return


if __name__ == "__main__":
    main()