#!/usr/bin/env python3 from setuptools import setup -name = 'amuse-fresco' -author = 'Steven Rieder and Inti Pelupessy' -author_email = '[email protected]' +name = "amuse-fresco" +author = "Steven Rieder and Inti Pelupessy" +author_email = "[email protected]" license_ = "Apache License 2.0" -url = 'https://github.com/rieder/fresco' +url = "https://github.com/rieder/fresco" install_requires = [ - 'wheel>=0.32', - 'amuse-framework>=2022.6.0', - 'scipy', - 'matplotlib', - 'astropy', + "wheel>=0.32", + "amuse-framework>=2022.6.0", + "scipy", + "matplotlib", + "astropy", ] setup_requires = [] -description = 'Make a realistic visualisation of a star cluster' +description = "Make a realistic visualisation of a star cluster" with open("README.md", "r", encoding="utf-8") as fh: long_description = fh.read() long_description_content_type = "text/markdown" extensions = [] all_data_files = [] -packages = [ - 'amuse.plot.fresco' -] +packages = ["amuse.plot.fresco"] -package_dir = { - 'amuse.plot.fresco': 'src/amuse/plot/fresco' -} +package_dir = {"amuse.plot.fresco": "src/amuse/plot/fresco"} -package_data = { -} +package_data = {} classifiers = [ # Maturity of amuse-fresco "Development Status :: 4 - Beta", # Intended audience

"Programming Language :: Python :: 3", "Programming Language :: Python :: 3.7", "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", - 'Topic :: Scientific/Engineering :: Astronomy', + "Topic :: Scientific/Engineering :: Astronomy", ] try: from src.amuse.plot.fresco.version import version + use_scm_version = False except ImportError: version = False - setup_requires += ['setuptools_scm',] + setup_requires += [ + "setuptools_scm", + ] use_scm_version = { "root": ".", "relative_to": __file__, "write_to": "src/amuse/plot/fresco/version.py", }

package_dir=package_dir, packages=packages, package_data=package_data, data_files=all_data_files, include_package_data=True, - python_requires='>=3.7, <4', + python_requires=">=3.7, <4", scripts=["bin/fresco.py"], )

# -*- coding: utf-8 -*- from __future__ import ( - print_function, - division, - ) + print_function, + division, +) import numpy as np from amuse.datamodel import Particles from amuse.units import units from amuse.ic.salpeter import new_salpeter_mass_distribution def new_field_stars( - N, - width=10 | units.parsec, - height=10 | units.parsec, - depth=100 | units.parsec, - massdistribution="salpeter", - agespread=3 | units.Gyr, - seed=1701, - ): + N, + width=10 | units.parsec, + height=10 | units.parsec, + depth=100 | units.parsec, + massdistribution="salpeter", + agespread=3 | units.Gyr, + seed=1701, +): np.random.seed(seed) stars = Particles(N) - stars.x = (np.random.random(N)-0.5) * width - stars.y = (np.random.random(N)-0.5) * height - stars.z = (np.random.random(N)-0.02) * depth + stars.x = (np.random.random(N) - 0.5) * width + stars.y = (np.random.random(N) - 0.5) * height + stars.z = (np.random.random(N) - 0.02) * depth if massdistribution == "salpeter": stars.mass = new_salpeter_mass_distribution(N) return stars

lambda x: B_lambda(x, T), ).in_(units.W / units.m**2) print(fb) print( - (fb * (1.0 | units.RSun) ** 2 / (1.0 | units.AU) ** 2).in_( - units.W / units.m**2 - ) + (fb * (1.0 | units.RSun) ** 2 / (1.0 | units.AU) ** 2).in_(units.W / units.m**2) ) print( ( energy_flux2(5778.0 | units.K) * (1.0 | units.RSun) ** 2

def calculate_effective_temperature(luminosity, radius): temp = ( np.nan_to_num( ( - (luminosity / (constants.four_pi_stefan_boltzmann * radius**2)) - ** 0.25 + (luminosity / (constants.four_pi_stefan_boltzmann * radius**2)) ** 0.25 ).value_in(units.K) ) | units.K ) return temp

import argparse import numpy as np import matplotlib -matplotlib.use('Agg') + +matplotlib.use("Agg") import matplotlib.pyplot as plt from amuse.units import units, nbody_system from amuse.datamodel import Particles from amuse.io import read_set_from_file

if parser is None: parser = argparse.ArgumentParser( formatter_class=argparse.ArgumentDefaultsHelpFormatter, ) parser.add_argument( - '--filetype', - dest='filetype', - default='amuse', - help='file type, valid are all types AMUSE can read', - ) - parser.add_argument( - '-s', - dest='starsfilename', - help='file containing stars (optional)', - ) - parser.add_argument( - '-g', - dest='gasfilename', - help='file containing gas (optional)', - ) - parser.add_argument( - '-f', - dest='followfilename', - help=( - 'file containing star keys to center on (optional, implies --com)' - ), - ) - parser.add_argument( - '-o', - dest='imagefilename', - help='write image to this file', - ) - parser.add_argument( - '--imagetype', - dest='imagetype', - default='png', - help='image file type', - ) - parser.add_argument( - '-b', - dest='sourcebands', - default='ubvri', - help='colour bands to use', - ) - parser.add_argument( - '-a', - dest='age', + "--filetype", + dest="filetype", + default="amuse", + help="file type, valid are all types AMUSE can read", + ) + parser.add_argument( + "-s", + dest="starsfilename", + help="file containing stars (optional)", + ) + parser.add_argument( + "-g", + dest="gasfilename", + help="file containing gas (optional)", + ) + parser.add_argument( + "-f", + dest="followfilename", + help=("file containing star keys to center on (optional, implies --com)"), + ) + parser.add_argument( + "-o", + dest="imagefilename", + help="write image to this file", + ) + parser.add_argument( + "--imagetype", + dest="imagetype", + default="png", + help="image file type", + ) + parser.add_argument( + "-b", + dest="sourcebands", + default="ubvri", + help="colour bands to use", + ) + parser.add_argument( + "-a", + dest="age", default=100 | units.Myr, type=units.Myr, - help='age of the stars', - ) - parser.add_argument( - '-w', - dest='width', + help="age of the stars", + ) + parser.add_argument( + "-w", + dest="width", default=5 | units.pc, type=units.pc, - help='image width', - ) - parser.add_argument( - '-x', - dest='plot_axes', - action='store_true', + help="image width", + ) + parser.add_argument( + "-x", + dest="plot_axes", + action="store_true", default=False, - help='plot axes', - ) - parser.add_argument( - '--ext', - dest='calculate_extinction', - action='store_true', + help="plot axes", + ) + parser.add_argument( + "--ext", + dest="calculate_extinction", + action="store_true", default=False, - help='include extinction by dust', - ) - parser.add_argument( - '--seed', - dest='seed', + help="include extinction by dust", + ) + parser.add_argument( + "--seed", + dest="seed", default=1701, type=int, - help='random seed', - ) - parser.add_argument( - '--vmax', - dest='vmax', + help="random seed", + ) + parser.add_argument( + "--vmax", + dest="vmax", default=0, type=float, - help='vmax value', - ) - parser.add_argument( - '--field', - dest='n_fieldstars', + help="vmax value", + ) + parser.add_argument( + "--field", + dest="n_fieldstars", default=0, type=int, - help='add N field stars (optional)', - ) - parser.add_argument( - '--ax', - dest='angle_x', + help="add N field stars (optional)", + ) + parser.add_argument( + "--ax", + dest="angle_x", default=0 | units.deg, type=units.deg, - help='Rotation step around x-axis', - ) - parser.add_argument( - '--ay', - dest='angle_y', + help="Rotation step around x-axis", + ) + parser.add_argument( + "--ay", + dest="angle_y", default=0 | units.deg, type=units.deg, - help='Rotation step around y-axis', - ) - parser.add_argument( - '--az', - dest='angle_z', + help="Rotation step around y-axis", + ) + parser.add_argument( + "--az", + dest="angle_z", default=0 | units.deg, type=units.deg, - help='Rotation step around z-axis', - ) - parser.add_argument( - '--frames', - dest='frames', + help="Rotation step around z-axis", + ) + parser.add_argument( + "--frames", + dest="frames", default=1, type=int, - help='Number of frames (>1: rotate around x,y,z)', - ) - parser.add_argument( - '--px', - dest='pixels', + help="Number of frames (>1: rotate around x,y,z)", + ) + parser.add_argument( + "--px", + dest="pixels", default=2048, type=int, - help='Number of pixels along each axis', - ) - parser.add_argument( - '--psf', - dest='psf_type', - default='hubble', + help="Number of pixels along each axis", + ) + parser.add_argument( + "--psf", + dest="psf_type", + default="hubble", help=( - 'PSF type. Looks for a .fits file of the given name, uses this if ' - 'it exists.\n' + "PSF type. Looks for a .fits file of the given name, uses this if " + "it exists.\n" 'Otherwise, "hubble", "wfc3", "wfpc2" and "gaussian" are valid ' - 'options.' - ) - ) - parser.add_argument( - '--sigma', - dest='psf_sigma', + "options." + ), + ) + parser.add_argument( + "--sigma", + dest="psf_sigma", default=1.0, type=float, - help='PSF sigma (if PSF type is gaussian)', - ) - parser.add_argument( - '--fl', - dest='fixed_luminosity', - action='store_true', + help="PSF sigma (if PSF type is gaussian)", + ) + parser.add_argument( + "--fl", + dest="fixed_luminosity", + action="store_true", default=False, - help='Use a fixed, equal luminosity and temperature for all stars', - ) - parser.add_argument( - '--contours', - dest='contours', - action='store_true', + help="Use a fixed, equal luminosity and temperature for all stars", + ) + parser.add_argument( + "--contours", + dest="contours", + action="store_true", default=False, - help='Plot gas contour lines', - ) - parser.add_argument( - '--com', - dest='use_com', - action='store_true', + help="Plot gas contour lines", + ) + parser.add_argument( + "--com", + dest="use_com", + action="store_true", default=False, - help='Center on center of mass', - ) - parser.add_argument( - '--xo', - dest='x_offset', + help="Center on center of mass", + ) + parser.add_argument( + "--xo", + dest="x_offset", default=0.0 | units.pc, type=units.pc, - help='X offset', - ) - parser.add_argument( - '--yo', - dest='y_offset', + help="X offset", + ) + parser.add_argument( + "--yo", + dest="y_offset", default=0.0 | units.pc, type=units.pc, - help='Y offset', - ) - parser.add_argument( - '--zo', - dest='z_offset', + help="Y offset", + ) + parser.add_argument( + "--zo", + dest="z_offset", default=0.0 | units.pc, type=units.pc, - help='Z offset', + help="Z offset", ) return parser.parse_args() def main():

# Derived settings image_size = [pixels, pixels] # If the nr of pixels is changed, zoom the PSF accordingly. - zoom_factor = pixels / 2048. + zoom_factor = pixels / 2048.0 if starsfilename: stars = read_set_from_file( starsfilename, filetype, close_file=True, ) if stellar_evolution and (age > 0 | units.Myr): - print(( - "Calculating luminosity/temperature for %s old stars..." - % (age) - )) + print(("Calculating luminosity/temperature for %s old stars..." % (age))) evolve_to_age(stars, age, stellar_evolution=se_code) elif args.fixed_luminosity: for band in sourcebands: setattr(stars, band + "_band", 1.0 | units.LSun) if use_com: if followfilename is not None: followstars = read_set_from_file( - followfilename, filetype, close_file=True, + followfilename, + filetype, + close_file=True, ) center_on_these_stars = followstars.get_intersecting_subset_in( stars, ) else:

if image_width == "max" | units.pc: minx = stars.x.min() maxx = stars.x.max() miny = stars.y.min() maxy = stars.y.max() - image_width = max(maxx-minx, maxy-miny) - x_offset = (maxx+minx) / 2 - y_offset = (maxy+miny) / 2 + image_width = max(maxx - minx, maxy - miny) + x_offset = (maxx + minx) / 2 + y_offset = (maxy + miny) / 2 stars.x -= x_offset stars.y -= y_offset stars.z -= z_offset - minx = (1.5 * -image_width/2) - maxx = (1.5 * image_width/2) - miny = (1.5 * -image_width/2) - maxy = (1.5 * image_width/2) + minx = 1.5 * -image_width / 2 + maxx = 1.5 * image_width / 2 + miny = 1.5 * -image_width / 2 + maxy = 1.5 * image_width / 2 stars = stars[stars.x > minx] stars = stars[stars.x < maxx] stars = stars[stars.y > miny] stars = stars[stars.y < maxy] else:

fieldstars = new_field_stars( n_fieldstars, width=image_width, height=image_width, ) - fieldstars.age = ( - minage - + ( - np.random.sample(n_fieldstars) - * (maxage - minage) - ) - ) + fieldstars.age = minage + (np.random.sample(n_fieldstars) * (maxage - minage)) evolve_to_age(fieldstars, 0 | units.yr, stellar_evolution=se_code) stars.add_particles(fieldstars) if gasfilename: gas = read_set_from_file(

# For gadget, need to divide by 2 to get the Fi value (??) # gas.h_smooth *= 0.5 # gas.radius = gas.h_smooth # Select only the relevant gas particles (plus a margin) - minx = (1.1 * -image_width/2) - maxx = (1.1 * image_width/2) - miny = (1.1 * -image_width/2) - maxy = (1.1 * image_width/2) + minx = 1.1 * -image_width / 2 + maxx = 1.1 * image_width / 2 + miny = 1.1 * -image_width / 2 + maxy = 1.1 * image_width / 2 gas_ = gas.select( - lambda x, y: - x > minx - and x < maxx - and y > miny - and y < maxy, - ["x", "y"] + lambda x, y: x > minx and x < maxx and y > miny and y < maxy, ["x", "y"] ) gas = gas_ else: gas = Particles() # gas.h_smooth = 0.05 | units.parsec

ax = fig.get_axes()[0] xmin, xmax = ax.get_xlim() ymin, ymax = ax.get_ylim() if not stars.is_empty(): - rotate(stars, (frames-1)*angle_x, (frames-1)*angle_y, (frames-1)*angle_z) + rotate( + stars, + (frames - 1) * angle_x, + (frames - 1) * angle_y, + (frames - 1) * angle_z, + ) if not gas.is_empty(): - rotate(gas, (frames-1)*angle_x, (frames-1)*angle_y, (frames-1)*angle_z) + rotate( + gas, (frames - 1) * angle_x, (frames - 1) * angle_y, (frames - 1) * angle_z + ) for frame in range(frames): frame += frames print( f"frame {frame}, " f"angle: {frame * angle_x.value_in(units.deg)} "

print(f"vmax = {vmax}") if not stars.is_empty(): ax.imshow( image, - origin='lower', + origin="lower", extent=[ xmin, xmax, ymin, ymax,

) gascontours[np.isnan(gascontours)] = 0.0 vmax = np.max(gascontours) / 2 # vmin = np.min(image[np.where(image > 0.0)]) vmin = vmax / 100 - levels = 10**( - np.linspace( - np.log10(vmin), - np.log10(vmax), - num=5, - ) - )[1:] + levels = ( + 10 + ** ( + np.linspace( + np.log10(vmin), + np.log10(vmax), + num=5, + ) + )[1:] + ) # print(vmin, vmax) # print(levels) ax.contour( - origin='lower', + origin="lower", levels=levels, colors="white", linewidths=0.1, extent=[ xmin,

image_size=image_size, ) ax.imshow( image, - origin='lower', + origin="lower", extent=[ xmin, xmax, ymin, ymax,