#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2023-02-13 16:20:51
"""Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution."""
import logging
import os
import fitsio
import numpy as np
from lsst_inaf_agile import util
[docs]
logger = logging.getLogger(__name__)
[docs]
class CatalogGalaxy:
"""Galaxy catalog class."""
def __init__(self, dirname: str, egg: dict):
"""Initialize the galaxy catalog."""
[docs]
self.bands = [s.strip() for s in self.egg["BANDS"][0]]
[docs]
self.rfbands = [s.strip() for s in self.egg["RFBANDS"][0]]
[docs]
self.catalog = self.get_catalog()
[docs]
def get_catalog(self, overwrite=False):
"""Generate the galaxy catalog and writes it to disk."""
# Create the Galaxy catalog
self.catalog = np.empty_like(self.egg["RA"][0], dtype=self.get_dtype())
filename = f"{self.dirname}/galaxy.fits"
if not os.path.exists(filename) or overwrite:
# column, type, description, unit
for c, t, _, _ in self.get_columns(self.bands, self.rfbands):
self.catalog[c] = self.get_galaxy(c).astype(t)
util.create_directory(filename)
fitsio.write(filename, self.catalog, clobber=True)
from lsst_inaf_agile.util import read_fits
self.catalog = read_fits(filename)
return self.catalog
@staticmethod
[docs]
def get_columns(bands: list[str], rfbands: list[str]) -> list[tuple]:
"""
Return the galaxy catalog columns, types, and descriptions.
The arguments 'bands' and 'rfbands' are used to generate the columns
for the fluxes and apparent magnitudes, respectively.
"""
ret = (
[
("ID", np.int64, "Unique ID", ""),
("RA", np.float64, "Right ascenscion", "deg"),
("DEC", np.float64, "Declination", "deg"),
("Z", np.float64, "Cosmological redshift", ""),
("D", np.float64, "Luminosity distance OR distance modulus for stars", "Mpc OR mag"),
("M", np.float64, "log10 of stellar mass", "Msun"),
("SFR", np.float64, "Star-formation rate", "Msun/yr"),
("PASSIVE", bool, "Is passive (non-star-forming)", ""),
("CLUSTERED", bool, "Is clustered", ""),
("BULGE_ANGLE", np.float64, "Galaxy bulge position angle", "deg"),
("BULGE_RADIUS", np.float64, "Galaxy bulge half-light radius", "arcsec"),
("BULGE_RATIO", np.float64, "Galaxy bulge ratio of minor to major axis", ""),
("DISK_ANGLE", np.float64, "Galaxy disk position angle", "deg"),
("DISK_RADIUS", np.float64, "Galaxy disk half-light radius", "arcsec"),
("DISK_RATIO", np.float64, "Galaxy disk ratio of minor to major axis", ""),
("AVLINES_BULGE", np.float64, "Emission line extinction in the bulge", "mag"),
("AVLINES_DISK", np.float64, "Emission line extinction in the disk", "mag"),
]
+ [
(b.strip() + f"_{comp}", np.float64, f"Galaxy {comp} {b} flux", "uJy")
for comp in ["bulge", "disk"]
for b in bands
]
+ [
("magabs_" + b.strip(), np.float64, f"Galaxy {b} absolute magnitude", "ABmag")
for b in rfbands
]
)
return ret
[docs]
def get_dtype(self) -> np.dtype:
"""Return the numpy dtype corresponding to the columns."""
dtype = []
for n, t, _, _ in self.get_columns(self.bands, self.rfbands):
dtype += [(n, t)]
return np.dtype(dtype)
[docs]
def get_galaxy(self, key: str):
"""Return the galaxy property corresponding to 'key'."""
logger.info(f"Getting galaxy attribute {key}")
_bands = [b.strip() for b in self.bands]
key_without_suffix = key.replace("_disk", "").replace("_bulge", "")
try:
if key_without_suffix in self.bands:
# 20250514
# idx = _bands.index(key_without_suffix)
idx = self.bands.index(key_without_suffix)
key_egg = "FLUX"
if "bulge" in key:
key_egg += "_BULGE"
if "disk" in key:
key_egg += "_DISK"
return self.egg[key_egg][0, :, idx]
elif "magabs_" in key:
idx = self.rfbands.index(key.replace("magabs_", ""))
return self.egg["RFMAG"][0, :, idx]
except ValueError:
# Fluxes not generated... return zerovalues
return np.full_like(self.egg["Z"][0], np.inf if "magabs_" in key else 0.0)
# Handles stellar mass
if key == "M":
# NOTE: Salpeter to Chabrier IMF conversion see e.g. Grylls+20
return self.egg["M"][0] - 0.24
# Handles star-formation rate
if key == "SFR":
# Where necessary to convert SFRs from the literature from Chabrier
# or Kroupa IMFs to the Salpeter IMF, we divide by constant factors
# of 0.63 (Chabrier) or 0.67 (Kroupa).
# SFR_salpeter = SFR_chabrier / 0.63
# ==> SFR_chabrier = 0.63 * SFR_salpeter
# (https://ned.ipac.caltech.edu/level5/March14/Madau/Madau3.html)
return self.egg["SFR"][0] * 0.63
# Handles any other egg column
return self.egg[key][0]
[docs]
def __getitem__(self, key: str):
"""Return the galaxy catalog column corresponding to 'key'."""
return self.get_galaxy(key)