Mock Milky Way stellar catalogues

Generate and analyze mock catalogues of stars in the Milky Way Galaxy.

Background/foreground stars are generated using the Galaxia set of models (Sharma et al., 2011). It is recommended you install and use the fork of Galaxia created for use by this library, though the original may also work.

Galactic stellar substructures, generated by external libraries, can be injected into the background galaxy to form a complete mock catalogue. Currently, globular clusters and stellar streams (from a few sources) are supported.

Supports the generation of synthetic photometry for the catalogues of stars using relevant isochrones, and a variety of slicing and filtering functions, which can be used to simulate membership selections using various photometric systems.

Installation

The mockgal package can be easily installed from this repository using pip.

pip install git+ssh://[email protected]/nmdickson/MockGal.git

or this repo can be cloned locally and installed with pip.

Usage

For more detailed documentation, see the documentation. (TODO sphinx/autodoc)

The mock catalogues of mockgal are typically made up of two components; the background galaxy and N stellar substructures. These (separate objects) come together to create the MockCatalogue object, which then handles all property getting, filtering, slicing and observing, for all components together.

import mockgal

In this example we will create a mock catalogue of stars around the globular cluster NGC362.

Creating the Mock Catalogue

Generate Background Galaxy (using Galaxia)

The background galaxy is generated using "Galaxia", as described above. You must ensure that Galaxia is installed and the executable is reachable within your path (alternatively you may also pass the path to the executable to galaxia_path below). mockgal will handle all calls to Galaxia internally.

In order to run Galaxia, you must ensure that the model data directory is reachable, by setting the GALAXIA_DATADIR environment variable. If this variable is not set, an attempt will be made to fall back to a ./GalaxiaData directory in your current working directory.

Running Galaxia is based on a set of model parameters. See the Galaxia documentation for details on each possible parameter. A number of good defaults are provided in parameters.py, but certain parameters must be given each time. Namely, the survey location and area, the base photometric system and the output file names should be provided. To ensure you don't run out of memory, it is probably also a good idea to provide an upper bounding magnitude limit. Any parameter not provided explicitly will use defaults mentioned above.

Remember that these stars can be sliced more precisely later, so err on the larger side for all limits. In particular, the magnitude limit cuts applied by Galaxia internally do not account for extinction, so you will want to use the slicing provided by mockgal to do so.

# Survey design

ll, bb = 305.89473384, -44.88931771  # Survey centre [degrees, galactic]
area = 50  # Survey (circular) area [square degrees]

# Output files

output_dir = "./NGC362/results"
output_fn = 'galaxy_around_ngc362'

# Photometric system

psys = 'CASTOR'  # You must ensure these isochrones exist, under this name
mag0, clr0 = "CASTOR_u", "CASTOR_uv-CASTOR_u"
mag0_upperlim = 40

# Parameter dict

params = {
    "outputFile": output_fn,
    "outputDir": output_dir,
    "seed": np.random.randint(100),
    "longitude": ll,
    "latitude": bb,
    "surveyArea": area,
    "photoSys": psys,
    "magcolorNames": ','.join((mag0, clr0)),
    "appMagLimits[1]": mag0_upperlim,
}

BackgroundGalaxy will support reading an already generated galaxia output file, given the path to it, if you do not want to rerun the model each time.

To generate a new background galaxy catalogue, the generate classmethod is used. If you want to include more than just the single base photometric system given above, you can also include them here by name, and they will be appended at the end.

bg = mockgal.BackgroundGalaxy.generate(params, extra_photosys=['LSST',])

This may take some time, depending on the parameters given, and will create an output ebf file at the given output location.

If this process fails, it may do so silently and catastrophically. Try adding the arg process_kw={'capture_output': False} to help debugging (this will break other stuff, only use to debug a Galaxia failure).

Load Stellar Substructure

Next we must get the model of NGC362 we wish to implant at the centre of our catalogue.

For this example, we will use a sampled model of the cluster created using GCfit. See the GCfit documentation for more details how to do this, and what the following means:

import gcfit
import astropy.coordinates as coord

# Taken from Dickson et al. 2023
th362 = dict(
    W0=5.459, M=0.278, rh=3.445, ra=1.410, g=1.454, delta=0.481,
    s2=0.000, F=2.517, a1=0.673, a2=0.780, a3=3.012, BHret=4.938, d=8.848
)
fm = gcfit.FittableModel(th362, observations=gcfit.Observations('NGC362'))

# Taken from Baumgardt catalogue
RV = 223.12 << (u.km / u.s)
μα, μδ = 6.694 << (u.mas / u.yr), -2.535 << (u.mas / u.yr)
RA, DEC = (6.023792, -72.081306) << u.deg

centre = coord.SkyCoord(ra=RA, dec=DEC, distance=fm.d,
                        pm_ra_cosdec=μα, pm_dec=μδ, radial_velocity=RV)

gcmodel = fm.sample(centre=centre)

There are multiple ways to create a substructure class, specific to the type and source of said structure. In this case, we use the GCfitGlobularCluster class. This class, and all other substructure classes, generates (from the model) the positions and masses of all the stars, and then uses the provided photometric system to interpolate magnitudes from the relevant isochrones.

The easiest way to do this is to reuse the photometric system object from the background galaxy, since they must match anyways.

ngc362 = mockgal.GCfitGlobularCluster('NGC362', model, bg.photosys)

For way to initialize different kinds of structures, see the documentation.

Technically mockgal can support multiple substructures, of any type, in one mock catalogue, if you want to for some reason.

Create Final Mock Catalogue

Creating the overall mock catalogue is as easy as providing the above objects.

mc = mockgal.MockCatalogue(bg, [ngc362,])

If you want to support observational uncertainties on the proper motions and magnitudes, you can also provide a few key observing quantities here as well.

import astropy.units as u
mc = mockgal.MockCatalogue(bg, [ngc362,],
                           exptime=550 << u.s, nexp=4, baseline=4 << u.yr)

Note that this functionality is somewhat unstable, and only implemented for the CASTOR photosystem currently.

Accessing any available property from the mock catalogue will fetch and process that property from all components. If a desired property doesn't exist in all components, getting it here will fail. Note that most property names come from the Galaxia model output, and might not be intuitive. If in doubt, check the available bg.attrs().

lon, lat, distance = mc.glon, mc.glat, mc.rad
pm_loncosb, pm_lat, vlos = mc.pm_l, mc.pm_b, mc.v_los

Most properties can simply be accessed normally. There also exists some helper functions for getting magnitudes and colours easily. Magnitudes can also typically be accessed from the object natively, using some correct formatting.

absmag, appmag = mc.CASTOR_u, mc.CASTOR_u_app
obsmag = mc.CASTOR_u_obs  # if available

Most of these attributes won't be visible to the object until runtime, so for a list of available properties like magnitudes, you must refer to the relevant photosystem (mc.photosys).

The select_component method will return a mask that can be used to select only the stars from a certain structure (or the background) from the returned properties.

Slicing and plotting

All stars in the mock catalogue can be sliced and filtered out, based on a number of criterion. See the documentation of MockCatalogue for details on all available slicing methods.

An initial slice will usually be applied to all catalogues at initialization, based on the magnitude slices in the Galaxia parameters, in order to correctly account for extinction, which Galaxia does not.

Slices are handled through the use of masked arrays on all fetched attributes and properties, and should be applied using one of the provided functions.

# Slice on some apparent magnitude depths

band_names = ('CASTOR_uv', 'CASTOR_u', 'CASTOR_g',
              'CASTOR_uv_split_bb', 'CASTOR_u_split_bb')
depths = (28.25, 28.00, 27.50, 27.40, 27.70)

mc.appmag_slice(band_names, upper_lim=depths)

# Slice out a small circle surrounding the cluster center

mc.field_slice(ngc362.get_bounding_circle(radius=2 * u.deg))

It's important to note here that apparent magnitude slices are handled slightly differently than most other slices, as you can slice (as done above) different bands to different depths. While most slicing will completely mask out the star in all properties, slicing on the apparent magnitudes, if the star is still visible in at least one band, will only mask the magnitude properties, and leave the star unmasked in all others.

mockgal supports a number of plots of the catalogue stars, with a variety of options to make slicing and membership selection tests easier.

Note that it may be necessary to really reduce the alpha value to see any structure, and if plotting all stars (i.e. not using the step arg) a lot of memory may be used.

import matplotlib.pyplot as plt

fig, (ax0, ax1) = plt.subplots(ncols=2, constrained_layout=True)

mc.plot_positions(fig=fig, ax=ax0, alpha=0.01, sep_comps=False, step=50)
mc.plot_positions(fig=fig, ax=ax1, alpha=0.01, sep_comps=True, step=50)

Often these plots can be used to provide useful visual hints about further slices necessary to perform something like membership selection. As an example, below, we use the property_2d_slice and CMD_slice methods to show how CASTOR may be used in the future to help with cluster membership selection.