CatalogBase.py

import numpy as np
import galsim
import copy
import cmath

from astropy.table import Table
from abc import abstractmethod, ABCMeta
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getABMAG


class CatalogBase(metaclass=ABCMeta):
    def __init__(self):
        pass

    @abstractmethod
    def load_sed(self, obj, **kward):
        pass

    @abstractmethod
    def _load(self, **kward):
        pass

    @abstractmethod
    def load_norm_filt(self, obj):
        return None
    
    @staticmethod
    def initialize_param():
        param = {
            "star":-1,
            "id":-1,
            "ra":0,
            "dec":0.,
            "ra_orig":0.,
            "dec_orig":0.,
            "z":0.,
            "sed_type":-1,
            "model_tag":"unknown",
            "mag_use_normal":100.,
            "theta":0.,
            "kappa":0.,
            "g1":0.,
            "g2":0.,
            "bfrac":0.,
            "av":0.,
            "redden":0.,
            "hlr_bulge":0.,
            "hlr_disk":0.,
            "ell_bulge":0.,
            "ell_disk":0.,
            "ell_tot":0.,
            "e1_disk":0.,
            "e2_disk":0.,
            "e1_bulge":0.,
            "e2_bulge":0.,
            "teff":0.,
            "logg":0.,
            "feh":0.,
            # C6 galaxies parameters
            "e1":0.,
            "e2":0.,
            "bulgemass":0.,
            "diskmass":0.,
            "size":0.,
            "detA":0.,
            "type":0,
            "veldisp":0.,
            "coeff": np.zeros(20),
            # Astrometry related
            "pmra":0.,
            "pmdec":0.,
            "rv":0.,
            "parallax":1e-9
        }
        return param

    @staticmethod
    def rotate_ellipticity(e1, e2, rotation=0., unit='radians'):
        if unit == 'degree':
            rotation = np.radians(rotation)
        e_total = np.sqrt(e1**2 + e2**2)
        phi = cmath.phase(complex(e1, e2))
        e1 = e_total * np.cos(phi + 2*rotation)
        e2 = e_total * np.sin(phi + 2*rotation)
        return e1, e2, e_total
    
    @staticmethod
    def convert_sed(mag, sed, target_filt, norm_filt=None):
        bandpass = target_filt.bandpass_full
        
        if norm_filt is not None:
            norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
        else:
            norm_filt = Table(
                np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
            )
            norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
        
        sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
                spectrum=sed,
                norm_thr=norm_filt,
                sWave=np.floor(norm_filt[norm_thr_rang_ids][0][0]),
                eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
        sed_photon = copy.copy(sed)
        sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
        sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
        # Get magnitude
        sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
        interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
        mag_csst = getABMAG(
            interFlux=interFlux,
            bandpass=bandpass
        )
        if target_filt.survey_type == "photometric":
            return sed_photon, mag_csst, interFlux
        elif target_filt.survey_type == "spectroscopic":
            del sed_photon
            return sed, mag_csst, interFlux