#!/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,