#!/usr/bin/env python3 from setuptools import setup -name = 'amuse-fresco' -author = 'Steven Rieder and Inti Pelupessy' -author_email = 'steven+fresco@rieder.nl' +name = "amuse-fresco" +author = "Steven Rieder and Inti Pelupessy" +author_email = "steven+fresco@rieder.nl" 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"], )

"""Functions for calculating a black body spectrum. """ from __future__ import ( - print_function, - division, - ) + print_function, + division, +) import numpy from numpy import ( - pi, e, + pi, + e, ) from amuse.units import units from amuse.units.constants import ( - kB, h, c, + kB, + h, + c, ) from amuse.units.constants import Rydberg_constant as Ry from amuse.units.constants import Stefan_hyphen_Boltzmann_constant as sigma def B_nu(nu, t): - return 2*h*nu**3/c**2 * 1. / (e**(h*nu/kB/t)-1) + return 2 * h * nu**3 / c**2 * 1.0 / (e ** (h * nu / kB / t) - 1) def B_lambda(wavelength, temperature): x = ( - 2 * h * c**2 + 2 + * h + * c**2 / wavelength**5 / ( - e**( - h * c + e + ** ( + h + * c / ( - (numpy.outer( - temperature.value_in(units.K), - wavelength.value_in(units.m) - ) | units.K * units.m) + ( + numpy.outer( + temperature.value_in(units.K), wavelength.value_in(units.m) + ) + | units.K * units.m + ) * kB - ) - ) - 1 + ) ) + - 1 ) + ) return x def energy(nu): - return h*nu + return h * nu def freq(e): - return e/h + return e / h def freq_from_wavenumber(k): - return c*k + return c * k def wavelength(nu): - return c/nu + return c / nu def freq_from_wavelength(wavelength): - return c/wavelength + return c / wavelength def wiens_lambda_max(T): b = 2897768.6 | units.nano(units.m) * units.K - return b/T + return b / T def wiens_T_from_lambda_max(wavelength): b = 2897768.6 | units.nano(units.m) * units.K - return b/wavelength + return b / wavelength -def energy_flux(T, lowfreq=0. | units.s**-1, N=1000): - nu = (numpy.arange(N+1)+1.)/N*(kB*T)/h*25. + lowfreq - b = pi*B_nu(nu, T) - return (b[1:]+b[:-1]).sum()/2*(nu[1]-nu[0]) +def energy_flux(T, lowfreq=0.0 | units.s**-1, N=1000): + nu = (numpy.arange(N + 1) + 1.0) / N * (kB * T) / h * 25.0 + lowfreq + b = pi * B_nu(nu, T) + return (b[1:] + b[:-1]).sum() / 2 * (nu[1] - nu[0]) -def energy_flux2(T, lambdas=None, throughput=1., N=1000): +def energy_flux2(T, lambdas=None, throughput=1.0, N=1000): if lambdas is None: lmax = wiens_lambda_max(T) - lambdas = lmax * 10**(-2. + 4. * numpy.arange(N+1)/float(N)) - b = pi*throughput*B_lambda(lambdas, T) - return numpy.trapz(b.number, x=lambdas.number) | (b.unit*lambdas.unit) + lambdas = lmax * 10 ** (-2.0 + 4.0 * numpy.arange(N + 1) / float(N)) + b = pi * throughput * B_lambda(lambdas, T) + return numpy.trapz(b.number, x=lambdas.number) | (b.unit * lambdas.unit) -def photon_flux(T, lowfreq=0. | units.s**-1, N=100000): - nu = (numpy.arange(N+1)+1.)/N*(kB*T)/h*25. + lowfreq - n = pi*B_nu(nu, T)/energy(nu) - return (n[1:]+n[:-1]).sum()/2*(nu[1]-nu[0]) +def photon_flux(T, lowfreq=0.0 | units.s**-1, N=100000): + nu = (numpy.arange(N + 1) + 1.0) / N * (kB * T) / h * 25.0 + lowfreq + n = pi * B_nu(nu, T) / energy(nu) + return (n[1:] + n[:-1]).sum() / 2 * (nu[1] - nu[0]) def total_bolometric_flux(T): - return sigma*T**4 + return sigma * T**4 if __name__ == "__main__": - T = 10000. | units.K + T = 10000.0 | units.K print((wiens_lambda_max(T))) print((energy_flux(T).in_(units.W * units.m**-2))) print((energy_flux2(T).in_(units.W * units.m**-2))) - print(((sigma*T**4).in_(units.W * units.m**-2))) + print(((sigma * T**4).in_(units.W * units.m**-2))) print() nf = photon_flux(T, lowfreq=freq_from_wavenumber(Ry)) print((numpy.log10(nf.value_in(units.cm**-2 * units.s**-1)))) print() a = photon_flux(T) print((numpy.log10(a.value_in(units.cm**-2 * units.s**-1)))) - b = sigma*T**4 / (kB*T)/2.7 + b = sigma * T**4 / (kB * T) / 2.7 print((numpy.log10(b.value_in(units.cm**-2 * units.s**-1)))) - print((b/a)) - print((nf/b)) + print((b / a)) + print((nf / b)) - print((wiens_T_from_lambda_max(300. | units.nano(units.m)))) - print((wiens_T_from_lambda_max(610. | units.nano(units.m)))) - print((wiens_T_from_lambda_max(920. | units.nano(units.m)))) + print((wiens_T_from_lambda_max(300.0 | units.nano(units.m)))) + print((wiens_T_from_lambda_max(610.0 | units.nano(units.m)))) + print((wiens_T_from_lambda_max(920.0 | units.nano(units.m))))

lambda x: x**3, lambda x: x**4, ] base2 = [ - lambda x: 1/(x+0.01)**3, - lambda x: x**0., - lambda x: 1/(1.01-x)**3, + lambda x: 1 / (x + 0.01) ** 3, + lambda x: x**0.0, + lambda x: 1 / (1.01 - x) ** 3, lambda x: x**3, lambda x: x**4, ] base3 = [ - lambda x: B_lambda(x, 5000. | units.K).number, - lambda x: B_lambda(x, 20000. | units.K).number, - lambda x: B_lambda(x, 2000. | units.K).number, - lambda x: B_lambda(x, 3000. | units.K).number, - lambda x: B_lambda(x, 10000. | units.K).number, + lambda x: B_lambda(x, 5000.0 | units.K).number, + lambda x: B_lambda(x, 20000.0 | units.K).number, + lambda x: B_lambda(x, 2000.0 | units.K).number, + lambda x: B_lambda(x, 3000.0 | units.K).number, + lambda x: B_lambda(x, 10000.0 | units.K).number, ] default_base = base1 class ColorConverter(object): def __init__( - self, source, target, - base=None, lmin=None, lmax=None, N=1000, + self, + source, + target, + base=None, + lmin=None, + lmax=None, + N=1000, ): if base is None: base = default_base if len(source) > len(base): raise Exception("provide enough base functions") if lmin is None: lmin = min( - [min(x['wavelength']) for x in source] - + [min(x['wavelength']) for x in target] + [min(x["wavelength"]) for x in source] + + [min(x["wavelength"]) for x in target] ) if lmax is None: lmax = max( - [max(x['wavelength']) for x in source] - + [max(x['wavelength']) for x in target] + [max(x["wavelength"]) for x in source] + + [max(x["wavelength"]) for x in target] ) self.lmin = lmin self.lmax = lmax self.dim = len(source) - self.base = base[0:self.dim] + self.base = base[0 : self.dim] self.source = source self.target = target self.dim2 = len(target) self.N = N - larray = lmin+(lmax-lmin)*numpy.array(list(range(N+1)))/N + larray = lmin + (lmax - lmin) * numpy.array(list(range(N + 1))) / N Amatrix = numpy.zeros((self.dim, self.dim)) for j, src in enumerate(source): - xp = src['wavelength'] - fp = src['throughput'] + xp = src["wavelength"] + fp = src["throughput"] f = numpy.interp( larray.number, xp=xp.value_in(larray.unit), fp=fp, - left=0., - right=0., + left=0.0, + right=0.0, ) for i, b in enumerate(self.base): - bint = f*self.fbase(b)(larray) + bint = f * self.fbase(b)(larray) Amatrix[j, i] = numpy.trapz(bint, x=larray.number) Bmatrix = numpy.zeros((self.dim2, self.dim)) for j, src in enumerate(target): - xp = src['wavelength'] - fp = src['throughput'] + xp = src["wavelength"] + fp = src["throughput"] f = numpy.interp( larray.number, xp=xp.value_in(larray.unit), fp=fp, - left=0., - right=0., + left=0.0, + right=0.0, ) for i, b in enumerate(self.base): - bint = f*self.fbase(b)(larray) + bint = f * self.fbase(b)(larray) Bmatrix[j, i] = numpy.trapz( bint, x=larray.number, )

self.Bmatrix = Bmatrix self.Ainv = numpy.linalg.inv(Amatrix) self.conversion_matrix = Bmatrix.dot(self.Ainv) def fbase(self, b): - return lambda x: b((x-self.lmin)/(self.lmax-self.lmin)) + return lambda x: b((x - self.lmin) / (self.lmax - self.lmin)) def convert(self, x): return self.conversion_matrix.dot(x)

class sRGB_linear_to_sRGB(object): def __init__(self): self.a = 0.055 def clip(self, x): - return x.clip(0., 1.) + return x.clip(0.0, 1.0) def convert(self, x): x = self.clip(x) a = numpy.where(x <= 0.0031308) b = numpy.where(x > 0.0031308) - x[a] = 12.92*x[a] - x[b] = (1+self.a) * x[b]**(1/2.4) - self.a + x[a] = 12.92 * x[a] + x[b] = (1 + self.a) * x[b] ** (1 / 2.4) - self.a return x

# -*- 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

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

"Parse command line arguments" 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 (amuse,starlab,txt,...)', - ) - parser.add_argument( - '-s', - dest='starsfilename', - default='', - help='file containing stars (optional)', - ) - parser.add_argument( - '-g', - dest='gasfilename', - default='', - help='file containing gas (optional)', - ) - parser.add_argument( - '-f', - dest='followfilename', + "--filetype", + dest="filetype", + default="amuse", + help="file type, valid are all types AMUSE can read (amuse,starlab,txt,...)", + ) + parser.add_argument( + "-s", + dest="starsfilename", + default="", + help="file containing stars (optional)", + ) + parser.add_argument( + "-g", + dest="gasfilename", + default="", + help="file containing gas (optional)", + ) + parser.add_argument( + "-f", + dest="followfilename", default=None, - help=\ - 'file containing star keys to center on (optional, implies --com)', - ) - parser.add_argument( - '-o', - dest='imagefilename', + help="file containing star keys to center on (optional, implies --com)", + ) + parser.add_argument( + "-o", + dest="imagefilename", default=None, - 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., - type=float, - help='age of the stars in Myr', - ) - parser.add_argument( - '-w', - dest='width', - default=5., - type=float, - help='image width in parsec', - ) - parser.add_argument( - '-x', - dest='plot_axes', - action='store_true', + 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.0, + type=float, + help="age of the stars in Myr", + ) + parser.add_argument( + "-w", + dest="width", + default=5.0, + type=float, + help="image width in parsec", + ) + 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, type=float, - help='Rotation step around x-axis in deg', - ) - parser.add_argument( - '--ay', - dest='angle_y', + help="Rotation step around x-axis in deg", + ) + parser.add_argument( + "--ay", + dest="angle_y", default=0, type=float, - help='Rotation step around y-axis in deg', - ) - parser.add_argument( - '--az', - dest='angle_z', + help="Rotation step around y-axis in deg", + ) + parser.add_argument( + "--az", + dest="angle_z", default=0, type=float, - help='Rotation step around z-axis in deg', - ) - parser.add_argument( - '--frames', - dest='frames', + help="Rotation step around z-axis in deg", + ) + 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', + ), + ) + 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, type=float, - help='X offset (in parsec)', - ) - parser.add_argument( - '--yo', - dest='y_offset', + help="X offset (in parsec)", + ) + parser.add_argument( + "--yo", + dest="y_offset", default=0.0, type=float, - help='Y offset (in parsec)', - ) - parser.add_argument( - '--zo', - dest='z_offset', + help="Y offset (in parsec)", + ) + parser.add_argument( + "--zo", + dest="z_offset", default=0.0, type=float, - help='Z offset (in parsec)', + help="Z offset (in parsec)", ) 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:

stars.z -= z_offset # Select only the relevant gas particles (plus a margin) # note: the margin should at least be half the PSF width - probably # more - 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

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,

"r": "bess-r.pass", "i": "bess-i.pass", } # Bessell 1990 filters -def get_filter_data( - instrument="WFPC_II_WFC3", - filters=bessellfilters -): +def get_filter_data(instrument="WFPC_II_WFC3", filters=bessellfilters): this_dir, this_filename = os.path.split(__file__) if this_dir == "": this_dir = "." data_dir = this_dir + "/data/" + instrument + "/"

) return filter_data def filter_band_flux( - fdata, - source, - wavelength=None, + fdata, + source, + wavelength=None, ): if fdata.__class__.__name__ == "str": fdata = get_filter_data()[fdata] if wavelength is None: - wavelength = fdata['wavelength'] - throughput = fdata['throughput'] + wavelength = fdata["wavelength"] + throughput = fdata["throughput"] else: - xp = fdata['wavelength'] - fp = fdata['throughput'] + xp = fdata["wavelength"] + fp = fdata["throughput"] throughput = numpy.interp( wavelength.value_in(units.angstrom), xp=xp.value_in(units.angstrom), fp=fp, - left=0., - right=0., + left=0.0, + right=0.0, ) src = throughput * source(wavelength) src = quantities.to_quantity(src) return numpy.trapz( src.number, x=wavelength.number, - ) | (src.unit*wavelength.unit) + ) | (src.unit * wavelength.unit) def filter_band_lambda(fdata): return ( - filter_band_flux(fdata, lambda x: x) - / filter_band_flux(fdata, lambda x: 1.) + filter_band_flux(fdata, lambda x: x) / filter_band_flux(fdata, lambda x: 1.0) ).in_(units.angstrom) def plot_filters( filters=bessellfilters,

ax = fig.add_subplot(1, 1, 1) filter_data = get_filter_data(filters=filters) n = 1000 for band in filters: - wavelength = ( - min_wavelength + (max_wavelength-min_wavelength) - * ( - numpy.array(range(n+1)) - / n - ) + wavelength = min_wavelength + (max_wavelength - min_wavelength) * ( + numpy.array(range(n + 1)) / n ) - xp = filter_data[filters[band]]['wavelength'] - fp = filter_data[filters[band]]['throughput'] + xp = filter_data[filters[band]]["wavelength"] + fp = filter_data[filters[band]]["throughput"] f = numpy.interp( wavelength.value_in(wavelength_unit), xp=xp.value_in(wavelength_unit), fp=fp, - left=0., - right=0., + left=0.0, + right=0.0, ) ax.plot(wavelength.value_in(wavelength_unit), f) - ax.set_xlabel(f'wavelength ({wavelength_unit})') - ax.set_ylabel('transmission') - + ax.set_xlabel(f"wavelength ({wavelength_unit})") + ax.set_ylabel("transmission") plt.show() if __name__ == "__main__": - plot_filters(filters=bessellfilters) - T = 5000. | units.K + T = 5000.0 | units.K fb = filter_band_flux( - filter_data["bess-u.pass"], - lambda x: B_lambda(x, T), - ).in_(units.W/units.m**2) + filter_data["bess-u.pass"], + lambda x: B_lambda(x, T), + ).in_(units.W / units.m**2) print(fb) print( - ( - fb * (1. | units.RSun)**2 - / (1. | 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. | units.K) - * (1. | units.RSun)**2 - / (1. | units.AU)**2 - ).in_(units.W/units.m**2) + energy_flux2(5778.0 | units.K) + * (1.0 | units.RSun) ** 2 + / (1.0 | units.AU) ** 2 + ).in_(units.W / units.m**2) )

""" import numpy as np from amuse.units import units -def map_to_grid( - *args, - **kwargs -): +def map_to_grid(*args, **kwargs): if len(args) == 2: return map_to_2d_grid(*args, **kwargs) return -1 def map_to_2d_grid( - x, - y, - weights=1, - image_size=(2048, 2048), - image_width=10 | units.parsec, - periodic=False, - mode="simple", + x, + y, + weights=1, + image_size=(2048, 2048), + image_width=10 | units.parsec, + periodic=False, + mode="simple", ): """ Returns a grid """ if mode == "simple":

) return -1 def map_to_2d_grid_simple( - x, - y, - weights=1, - image_size=(2048, 2048), - image_width=10 | units.parsec, + x, + y, + weights=1, + image_size=(2048, 2048), + image_width=10 | units.parsec, ): try: x_size = image_size[0] y_size = image_size[1] except TypeError:

x_0 = np.floor(x_px) y_0 = np.floor(y_px) x_1 = x_0 + 1 y_1 = y_0 + 1 - weight_x0 = (x_1 - x_px) - weight_y0 = (y_1 - y_px) - weight_x1 = (x_px - x_0) - weight_y1 = (y_px - y_0) + weight_x0 = x_1 - x_px + weight_y0 = y_1 - y_px + weight_x1 = x_px - x_0 + weight_y1 = y_px - y_0 pos_weights = [0, 0, 0, 0] pos_weights[0] = weight_x0 * weight_y0 pos_weights[1] = weight_x0 * weight_y1 pos_weights[2] = weight_x1 * weight_y0 pos_weights[3] = weight_x1 * weight_y1

def evolve_to_age(stars, age, stellar_evolution="SeBa"): "Evolve stars to specified age with specified code" if stellar_evolution == "SeBa": from amuse.community.seba.interface import SeBa + stellar_evolution = SeBa() elif stellar_evolution == "SSE": from amuse.community.sse.interface import SSE + stellar_evolution = SSE() # SSE can result in nan values for luminosity/radius else: raise ( f"No such stellar evolution code {stellar_evolution} or no code "

) stellar_evolution.particles.add_particles(stars) if age > 0 | units.yr: stellar_evolution.evolve_model(age) - stars.luminosity = np.nan_to_num( - stellar_evolution.particles.luminosity.value_in(units.LSun) - ) | units.LSun + stars.luminosity = ( + np.nan_to_num(stellar_evolution.particles.luminosity.value_in(units.LSun)) + | units.LSun + ) stars.radius = stellar_evolution.particles.radius # prevent zero/nan radius. location_of_nan_radius = np.where( - np.nan_to_num( - stars.radius.value_in(units.RSun) - ) == 0. + np.nan_to_num(stars.radius.value_in(units.RSun)) == 0.0 ) stars[location_of_nan_radius].radius = 0.01 | units.RSun stellar_evolution.stop() return def calculate_effective_temperature(luminosity, radius): - temp = np.nan_to_num( - ( + temp = ( + np.nan_to_num( ( - luminosity - / ( - constants.four_pi_stefan_boltzmann - * radius**2 - ) - )**.25 - ).value_in(units.K) - ) | units.K + (luminosity / (constants.four_pi_stefan_boltzmann * radius**2)) + ** 0.25 + ).value_in(units.K) + ) + | units.K + ) return temp def make_image( maps=None, stars=None, gas=None, converter=None, image_width=[ - 10. | units.parsec, - 10. | units.parsec, + 10.0 | units.parsec, + 10.0 | units.parsec, ], image_size=[1024, 1024], percentile=0.9995, - age=0. | units.Myr, + age=0.0 | units.Myr, sourcebands="ubvri", vmax=None, calc_temperature=True, mapper_code=None, # "FiMap" zoom_factor=1.0,

if extinction: # Extinction can currently only be handled with FiMap mapper_code = "FiMap" if mapper_code == "FiMap": + def mapper(): from amuse.community.fi.interface import FiMap + mapper = FiMap(converter, mode="openmp") # mapper.parameters.minimum_distance = 1. | units.AU mapper.parameters.image_size = image_size # mapper.parameters.image_target = image_target

# mapper.parameters.projection_mode = projection # mapper.parameters.image_angle = horizontal_angle # mapper.parameters.viewpoint = viewpoint mapper.parameters.extinction_flag = extinction return mapper + else: # Gridify as default mapper = None mapper_code = "gridify"

) return image def column_density_map( - gas, - image_width=10. | units.parsec, - image_size=[1024, 1024], - mapper_factory=None, - mapper_code=None, - zoom_factor=1.0, - psf_type="gaussian", - psf_sigma=10.0, - return_vmax=False, + gas, + image_width=10.0 | units.parsec, + image_size=[1024, 1024], + mapper_factory=None, + mapper_code=None, + zoom_factor=1.0, + psf_type="gaussian", + psf_sigma=10.0, + return_vmax=False, ): if mapper_code == "FiMap": if callable(mapper_factory): mapper = mapper_factory()

return im def image_from_stars( stars, - image_width=10. | units.parsec, + image_width=10.0 | units.parsec, image_size=[1024, 1024], percentile=0.9995, calc_temperature=True, - age=0. | units.Myr, + age=0.0 | units.Myr, sourcebands="ubvri", gas=None, vmax=None, mapper_factory=None, mapper_code=None,

psf_type=psf_type, psf_file=psf_file, psf_sigma=psf_sigma, ) if return_vmax: - return rgb['pixels'], vmax - return rgb['pixels'] + return rgb["pixels"], vmax + return rgb["pixels"] def initialise_image( - fig=None, - dpi=150, - image_size=[2048, 2048], - length_unit=units.parsec, - image_width=5 | units.parsec, - plot_axes=True, - subplot=0, - x_offset=0 | units.parsec, - y_offset=0 | units.parsec, - z_offset=0 | units.parsec, + fig=None, + dpi=150, + image_size=[2048, 2048], + length_unit=units.parsec, + image_width=5 | units.parsec, + plot_axes=True, + subplot=0, + x_offset=0 | units.parsec, + y_offset=0 | units.parsec, + z_offset=0 | units.parsec, ): if fig is None: if plot_axes: left = 0.2 bottom = 0.2 else: - left = 0. - bottom = 0. + left = 0.0 + bottom = 0.0 right = 1.0 top = 1.0 figwidth = image_size[0] / dpi / (right - left) figheight = image_size[1] / dpi / (top - bottom) figsize = (figwidth, figheight)

ax = fig.get_axes()[subplot] ax.cla() ax.set_xlabel("X (%s)" % (length_unit)) ax.set_ylabel("Y (%s)" % (length_unit)) ax.set_aspect(1) - ax.spines['top'].set_visible(False) - ax.spines['right'].set_visible(False) - ax.spines['bottom'].set_visible(False) - ax.spines['left'].set_visible(False) - ax.set_facecolor('black') + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.spines["bottom"].set_visible(False) + ax.spines["left"].set_visible(False) + ax.set_facecolor("black") return fig

x = numpy.array(x) y = numpy.array(y) z = numpy.array(z) xyz_data = dict() -xyz_data['x'] = dict(wavelength=cielam, throughput=x) -xyz_data['y'] = dict(wavelength=cielam, throughput=y) -xyz_data['z'] = dict(wavelength=cielam, throughput=z) +xyz_data["x"] = dict(wavelength=cielam, throughput=x) +xyz_data["y"] = dict(wavelength=cielam, throughput=y) +xyz_data["z"] = dict(wavelength=cielam, throughput=z)

# import logging # logging.basicConfig(level=logging.DEBUG) from .filters import ( - filter_band_flux, get_filter_data, filter_band_lambda, + filter_band_flux, + get_filter_data, + filter_band_lambda, ) from .xyz import xyz_data from .blackbody import B_lambda -from .color_converter import ( - ColorConverter, - XYZ_to_sRGB_linear, sRGB_linear_to_sRGB -) +from .color_converter import ColorConverter, XYZ_to_sRGB_linear, sRGB_linear_to_sRGB def Convolve( - image, - kernel, -): - result = convolve_fft(image, kernel, boundary='fill') + image, + kernel, +): + result = convolve_fft(image, kernel, boundary="fill") return result def get_psf( instrument="WFPC_II_WFC3",

psf_filename = datadir + band + "%2.2i.fits" % i f = pyfits.open(psf_filename) if zoom_factor != 1.0: psf[band + str(i)] = zoom( numpy.array(f[0].data), - zoom_factor, order=3, mode='constant', cval=0.0, + zoom_factor, + order=3, + mode="constant", + cval=0.0, prefilter=True, ) else: psf[band + str(i)] = numpy.array(f[0].data) return psf

mapper_gas, stars, gas, dust_to_gas=0.01, dust_extinction_cross_section_v_band=4.9e-22 | units.cm**2, - dust_albedo_v_band=.01, + dust_albedo_v_band=0.01, Nstar=50, # number of stars to use for reflected light on gas/dust ): - mapper_stars.weight = getattr(stars, band + "_band").value_in(units.LSun) if gas.is_empty(): return lambda_eff = filter_band_lambda("bess-" + band + ".pass") lambda_v = filter_band_lambda("bess-v.pass") f_H = 0.7 dust_cross_section = ( - dust_extinction_cross_section_v_band - * (lambda_eff / lambda_v)**-1.5 + dust_extinction_cross_section_v_band * (lambda_eff / lambda_v) ** -1.5 ) - dust_albedo = dust_albedo_v_band * (lambda_eff / lambda_v)**0. - - mapper_gas.opacity_area = ( - f_H - * gas.mass.value_in(units.amu) - * dust_cross_section - ) + dust_albedo = dust_albedo_v_band * (lambda_eff / lambda_v) ** 0.0 + + mapper_gas.opacity_area = f_H * gas.mass.value_in(units.amu) * dust_cross_section weight = numpy.zeros(len(mapper_gas)) stars_ordered_by_lum = stars.sorted_by_attribute(band + "_band") Nstar = min(Nstar, len(stars)) for star in stars_ordered_by_lum[-Nstar:]: d2 = ( - (gas.x - star.x)**2 - + (gas.y - star.y)**2 - + (gas.z - star.z)**2 + (gas.x - star.x) ** 2 + + (gas.y - star.y) ** 2 + + (gas.z - star.z) ** 2 + 0.25 * gas.radius**2 ) flux = getattr(star, band + "_band") / (4 * numpy.pi * d2) - weight += ( - flux - * mapper_gas.opacity_area - * dust_albedo - ).value_in(units.LSun) + weight += (flux * mapper_gas.opacity_area * dust_albedo).value_in(units.LSun) mapper_gas.weight = weight return def rgb_frame(

dryrun=False, vmax=None, percentile=0.9995, multi_psf=False, sourcebands="ubvri", - image_width=12. | units.parsec, + image_width=12.0 | units.parsec, image_size=[1024, 1024], mapper_factory=None, gas=None, mapper_code=None, zoom_factor=1.0, psf_type="hubble", psf_file=None, psf_sigma=1.0, verbose=False, ): - if gas is None: gas = Particles() if verbose: print("luminosities..") for band in sourcebands: - if not hasattr( - stars, - band + "_band" - ): + if not hasattr(stars, band + "_band"): setattr( stars, band + "_band", - 4 * numpy.pi * stars.radius**2 * - filter_band_flux( + 4 + * numpy.pi + * stars.radius**2 + * filter_band_flux( "bess-" + band + ".pass", lambda x: B_lambda(x, stars.temperature), ), ) if verbose:

mapper.stop() else: # Use simpler python mapping script from .gridify import map_to_grid + stars_in_mapper = stars.copy() gas_in_mapper = gas.copy() raw_images = dict() for band in sourcebands:

psf = get_psf(zoom_factor=zoom_factor) if multi_psf: a = numpy.arange(image_size[0]) / float(image_size[0] - 1) b = numpy.arange(image_size[1]) / float(image_size[1] - 1) w1 = numpy.outer(a, b) - w2 = numpy.outer(1. - a, b) - w3 = numpy.outer(a, 1. - b) - w4 = numpy.outer(1. - a, 1. - b) + w2 = numpy.outer(1.0 - a, b) + w3 = numpy.outer(a, 1.0 - b) + w4 = numpy.outer(1.0 - a, 1.0 - b) for key, val in list(raw_images.items()): - im1 = Convolve(val, psf[key + '0']) - im2 = Convolve(val, psf[key + '1']) - im3 = Convolve(val, psf[key + '2']) - im4 = Convolve(val, psf[key + '3']) - convolved_images[key] = ( - w1 * im1 - + w2 * im2 - + w3 * im3 - + w4 * im4 - ) + im1 = Convolve(val, psf[key + "0"]) + im2 = Convolve(val, psf[key + "1"]) + im3 = Convolve(val, psf[key + "2"]) + im4 = Convolve(val, psf[key + "3"]) + convolved_images[key] = w1 * im1 + w2 * im2 + w3 * im3 + w4 * im4 else: for key, val in list(raw_images.items()): - im1 = Convolve(val, psf[key + '0']) + im1 = Convolve(val, psf[key + "0"]) convolved_images[key] = im1 elif psf_type == "gaussian": for key, val in list(raw_images.items()): im1 = gaussian_filter(val, sigma=psf_sigma, order=0) convolved_images[key] = im1 if verbose: print("..conversion to rgb") filter_data = get_filter_data() - source = [filter_data['bess-' + x + '.pass'] for x in sourcebands] - - target = [xyz_data['x'], xyz_data['y'], xyz_data['z']] + source = [filter_data["bess-" + x + ".pass"] for x in sourcebands] + + target = [xyz_data["x"], xyz_data["y"], xyz_data["z"]] conv = ColorConverter(source, target) ubv = numpy.array([convolved_images[x] for x in sourcebands])