Catalog_example.py

import os
import numpy as np
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate

from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar


class Catalog(CatalogBase):
    """An user customizable class for reading in catalog(s) of objects and SEDs.

    NOTE: must inherit the "CatalogBase" abstract class

    ...

    Attributes
    ----------
    cat_dir : str
        a directory that contains the catalog file(s)
    star_path : str
        path to the star catalog file
    star_SED_path : str
        path to the star SED data
    objs : list
        a list of observation_sim.mock_objects (Star, Galaxy, or Quasar)
        NOTE: must have "obj" list when implement your own Catalog

    Methods
    ----------
    load_sed(obj, **kwargs):
        load the corresponding SED data for one object
    load_norm_filt(obj):
        load the filter throughput for the input catalog's photometric system.
    """

    def __init__(self, config, chip, **kwargs):
        """Constructor method.

        Parameters
        ----------
        config : dict
            configuration dictionary which is parsed from the input YAML file
        chip: observation_sim.instruments.Chip
            an observation_sim.instruments.Chip instance, can be used to identify the band etc.
        **kwargs : dict
            any other needed input parameters (in key-value pairs), please modify corresponding
            initialization call in "ObservationSim.py" as you need.

        Returns
        ----------
        None
        """

        super().__init__()
        self.cat_dir = os.path.join(
            config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
        self.chip = chip
        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"]:
            star_file = config["catalog_options"]["input_path"]["star_cat"]
            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
            self.star_path = os.path.join(self.cat_dir, star_file)
            self.star_SED_path = os.path.join(
                config["data_dir"], star_SED_file)
        # NOTE: must call _load() method here to read in all objects
        self.objs = []
        self._load()

    def _load(self, **kwargs):
        """Read in all objects in from the catalog file(s).

        This is a must implemented method which is used to read in all objects, and
        then convert them to observation_sim.mock_objects (Star, Galaxy, or Quasar).

        Currently, 
        the model of observation_sim.mock_objects.Star class requires:
            param["star"] : int
                specify the object type: 0: galaxy, 1: star, 2: quasar
            param["id"] : int
                ID number of the object
            param["ra"] : float
                Right ascension (in degrees)
            param["dec"] : float
                Declination (in degrees)
            param["mag_use_normal"]: float
                the absolute magnitude in a particular filter
                NOTE: if that filter is not the corresponding CSST filter, the 
                load_norm_filt(obj) function must be implemented to load the filter
                throughput of that particular photometric system

        the model of observation_sim.mock_objects.Galaxy class requires:
            param["star"] : int
                specify the object type: 0: galaxy, 1: star, 2: quasar
            param["id"] : int
                ID number of the object
            param["ra"] : float
                Right ascension (in degrees)
            param["dec"] : float
                Declination (in degrees)
            param["mag_use_normal"]: float
                the absolute magnitude in a particular filter
                NOTE: if that filter is not the corresponding CSST filter, the 
                load_norm_filt(obj) function must be implemented to load the filter
                throughput of that particular photometric system
            param["bfrac"] : float
                the bulge fraction
            param["hlr_bulge"] : float
                the half-light-radius of the bulge
            param["hlr_disk"] : float
                the half-light-radius of the disk
            param["e1_bulge"], param["e2_bulge"] : float
                the ellipticity of the bulge components
            param["e1_disk"], param["e2_disk"] : float
                the ellipticity of the disk components
            (Optional parameters):
            param['disk_sersic_idx']: float
                Sersic index for galaxy disk component
            param['bulge_sersic_idx']: float
                Sersic index for galaxy bulge component
            param['g1'], param['g2']: float
                Reduced weak lensing shear components (valid for shear type: catalog)
        the model of observation_sim.mock_objects.Galaxy class requires:
            Currently a Quasar is modeled as a point source, just like a Star.

        NOTE: To construct an object, according to its type, just call:
            Star(param), Galaxy(param), or Quasar(param)

        NOTE: All constructed objects should be appened to "self.objs".

        NOTE: Any other parameters can also be set within "param" dict:
            Used to calculate required quantities and/or SEDs etc.

        Parameters
        ----------
        **kwargs : dict
            any other needed input parameters (in key-value pairs), please modify corresponding
            initialization call in "ObservationSim.py" as you need.

        Returns
        ----------
        None
        """

        stars = Table.read(self.star_path)
        nstars = stars['sourceID'].size
        for istars in range(nstars):
            param = self.initialize_param()
            param['id'] = istars + 1
            param['ra'] = stars['RA'][istars]
            param['dec'] = stars['Dec'][istars]
            param['sed_type'] = stars['sourceID'][istars]
            param['model_tag'] = stars['model_tag'][istars]
            param['z'] = 0.0
            param['star'] = 1   # Star
            param['mag_use_normal'] = stars['app_sdss_g'][istars]
            obj = Star(param)
            self.objs.append(obj)

    def load_sed(self, obj, **kwargs):
        """Load the corresponding SED data for a particular obj.

        Parameters
        ----------
        obj : instance of a particular class in observation_sim.mock_objects
            the object to get SED data for
        **kwargs : dict
            other needed input parameters (in key-value pairs), please modify corresponding
            initialization call in "ObservationSim.py" as you need.

        Returns
        ----------
        sed : Astropy.Table
            the SED Table with two columns (namely, "WAVELENGTH", "FLUX"):
                sed["WAVELENGTH"] : wavelength in Angstroms
                sed["FLUX"] : fluxes in photons/s/m^2/A
            NOTE: the range of wavelengthes must at least cover [2450 - 11000] Angstorms
        """

        if obj.type == 'star':
            wave = Table.read(self.star_SED_path,
                              path=f"/SED/wave_{obj.model_tag}")
            flux = Table.read(self.star_SED_path, path=f"/SED/{obj.sed_type}")
            wave, flux = wave['col0'].data, flux['col0'].data
        else:
            raise ValueError("Object type not known")
        speci = interpolate.interp1d(wave, flux)
        lamb = np.arange(2400, 11001 + 0.5, 0.5)
        y = speci(lamb)
        # erg/s/cm^2/A --> photons/s/m^2/A
        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
        return sed

    def load_norm_filt(self, obj):
        """Load the corresponding thourghput for the input magnitude "param["mag_use_normal"]".

        NOTE: if the input magnitude is already in CSST magnitude, simply return None

        Parameters
        ----------
        obj : instance of a particular class in observation_sim.mock_objects
            the object to get thourghput data for

        Returns
        ----------
        norm_filt : Astropy.Table
            the throughput Table with two columns (namely, "WAVELENGTH", "SENSITIVITY"):
                norm_filt["WAVELENGTH"] : wavelengthes in Angstroms
                norm_filt["SENSITIVITY"] : efficiencies
        """
        return None