Galaxy.py

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

from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.MockObject import MockObject

class Galaxy(MockObject):
    def __init__(self, param, rotation=None):
        super().__init__(param)
        self.thetaR = self.param["theta"]
        self.bfrac = self.param["bfrac"]
        self.hlr_disk = self.param["hlr_disk"]
        self.hlr_bulge = self.param["hlr_bulge"]

        # Extract ellipticity components
        self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees)
        self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees)
        self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees)
        self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2
        self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2
        self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2

        if rotation is not None:
            self.rotateEllipticity(rotation)

    # def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
    #     if survey_type == "photometric":
    #         norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
    #         if sed_templates is None:
    #             # Read SED data directly
    #             itype = objtypes[cosids==self.sed_type][0]
    #             sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type))
    #             if not os.path.exists(sed_file):
    #                 raise ValueError("!!! No SED found.")
    #             sed_data = Table.read(sed_file, format="ascii")
    #             wave, flux = sed_data["observedLambda"].data, sed_data["observedFlux"].data
    #         else:
    #             # Load SED from templates
    #             sed_data = sed_templates[self.sed_type]
    #             # redshift, intrinsic extinction
    #             sed_data = getObservedSED(
    #                 sedCat=sed_data, 
    #                 redshift=self.z, 
    #                 av=self.param['av'], 
    #                 redden=self.param['redden'])
    #             wave, flux = sed_data[0], sed_data[1] 
    #         flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
    #         sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
    #         # Get scaling factor for SED
    #         sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
    #             spectrum=sed_photon,
    #             norm_thr=normFilter,
    #             sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
    #             eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
    #         sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
    #         # Convert to galsim.SED object
    #         spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
    #         self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
    #         # Get magnitude
    #         interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
    #         self.param['mag_%s'%target_filt.filter_type] = getABMAG(
    #             interFlux=interFlux, 
    #             bandpass=target_filt.bandpass_full)
    #         # print('mag_use_normal = ', self.param['mag_use_normal'])
    #         # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
    #         # print('redshift = %.3f'%(self.z))
    #         # print('sed_type = %d, av = %.2f, redden = %d'%(self.sed_type, self.param['av'], self.param['redden']))

    #     elif survey_type == "spectroscopic":
    #         if sed_templates is None:
    #             self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes)
    #         else:
    #             sed_data = sed_templates[self.sed_type]
    #             sed_data = getObservedSED(
    #                 sedCat=sed_data, 
    #                 redshift=self.z, 
    #                 av=self.param['av'], 
    #                 redden=self.param['redden'])
    #             speci = interpolate.interp1d(sed_data[0], sed_data[1])
    #             lamb = np.arange(2500, 10001 + 0.5, 0.5)
    #             y = speci(lamb)
    #             # erg/s/cm2/A --> photo/s/m2/A
    #             all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
    #             self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))


    def unload_SED(self):
        """(Test) free up SED memory
        """
        del self.sed

    # def sedPhotons(self, sed_path, cosids, objtypes):

    #     itype = objtypes[cosids == self.sed_type][0]
    #     sed_file = os.path.join(sed_path, itype + "_ID%s.sed" % (self.sed_type))
    #     if not os.path.exists(sed_file):
    #         raise ValueError("!!! No SED found.")
    #     sed = Table.read(sed_file, format="ascii")
    #     spec_data = {}
    #     f_orig = sed["observedFlux"].data
    #     w_orig = sed["observedLambda"].data

    #     speci = interpolate.interp1d(w_orig, f_orig)
    #     lamb = np.arange(2500, 10001 + 0.5, 0.5)
    #     y = speci(lamb)
    #     # erg/s/cm2/A --> photo/s/m2/A
    #     all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
    #     self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))

    def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
        if len(psf_list) != len(bandpass_list):
            raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
        objs = []
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
            return -1
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
                return -1
            
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                return -1

            psf = psf_list[i]
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
            gal = gal.withFlux(nphotons)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
            gal = galsim.Convolve(psf, gal)
            objs.append(gal)
        final = galsim.Sum(objs)
        return final

    def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150.):
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
            return False

        nphotons_sum = 0
        photons_list = []
        xmax, ymax = 0, 0

        # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge))
        big_galaxy = False
        if self.hlr_disk > 3.0: # Very big galaxy
            big_galaxy = True

        # (TEST) Galsim Parameters
        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)

        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]

            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
                # return False
                continue
            
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                # return False
                continue
            nphotons_sum += nphotons
            # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))

            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
            gal = gal.withFlux(nphotons)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)

            if self.hlr_disk < 10.0: # Not apply PSF for very big galaxy
                gal = galsim.Convolve(psf, gal)

            # Use (explicit) stamps to draw
            stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True)
            xmax = max(xmax, stamp.xmax)
            ymax = max(ymax, stamp.ymax)
            photons = stamp.photons
            photons.x += self.x_nominal
            photons.y += self.y_nominal
            photons_list.append(photons)

        # print('xmax = %d, ymax = %d '%(xmax, ymax))

        stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1))
        stamp.wcs = self.localWCS
        stamp.setCenter(self.x_nominal, self.y_nominal)
        bounds = stamp.bounds & chip.img.bounds
        stamp[bounds] = chip.img[bounds]

        if not big_galaxy:
            for i in range(len(photons_list)):
                if i == 0:
                    chip.sensor.accumulate(photons_list[i], stamp)
                else:
                    chip.sensor.accumulate(photons_list[i], stamp, resume=True)
        else:
            sensor = galsim.Sensor()
            for i in range(len(photons_list)):
                if i == 0:
                    sensor.accumulate(photons_list[i], stamp)
                else:
                    sensor.accumulate(photons_list[i], stamp, resume=True)

        # print(stamp.array.sum())
        # chip.img[bounds] += stamp[bounds]
        chip.img[bounds] = stamp[bounds]
        # print("nphotons_sum = ", nphotons_sum)
        del photons_list
        del stamp
        return True, pos_shear

    def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
                         exptime=150., normFilter=None, grating_split_pos=3685):

        norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
        sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
                                                      norm_thr=normFilter,
                                                      sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
                                                      eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
        if sedNormFactor == 0:
            return False
        normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
                          names=('WAVELENGTH', 'FLUX'))

        big_galaxy = False
        if self.hlr_disk > 3.0: # Very big galaxy
            big_galaxy = True

        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)
        # nphotons_sum = 0
        xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
        grating_split_pos_chip = chip.bound.xmin + grating_split_pos
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]

            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
            gal = gal.withFlux(tel.pupil_area * exptime)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
            gal = galsim.Convolve(psf, gal)

            starImg = gal.drawImage(wcs=self.localWCS)

            origin_star = [self.y_nominal - (starImg.center.y - starImg.ymin),
                           self.x_nominal - (starImg.center.x - starImg.xmin)]

            gal_origin = [origin_star[0], origin_star[1]]
            gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]

            if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
                subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
                ## part img disperse

                subImg_p1 = starImg.array[:, 0:subSlitPos]
                star_p1 = galsim.Image(subImg_p1)
                origin_p1 = origin_star
                xcenter_p1 = min(self.x_nominal,grating_split_pos_chip-1) - chip.bound.xmin
                ycenter_p1 = self.y_nominal-chip.bound.ymin

                sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
                                    ycenter=ycenter_p1, origin=origin_p1,
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[0],
                                    isAlongY=0)

                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus)

                subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
                star_p2 = galsim.Image(subImg_p2)
                origin_p2 = [origin_star[0], grating_split_pos_chip]
                xcenter_p2 = max(self.x_nominal, grating_split_pos_chip - 1) - chip.bound.xmin
                ycenter_p2 = self.y_nominal - chip.bound.ymin

                sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
                                       ycenter=ycenter_p2, origin=origin_p2,
                                       tar_spec=normalSED,
                                       band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                       conf=chip.sls_conf[1],
                                       isAlongY=0)

                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus)

                del sdp_p1
                del sdp_p2
            elif grating_split_pos_chip<=gal_origin[1]:
                sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin,
                                    ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star,
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[1],
                                    isAlongY=0)
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus)
                del sdp
            elif grating_split_pos_chip>=gal_end[1]:
                sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin,
                                    ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star,
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[0],
                                    isAlongY=0)
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus)
                del sdp

            # print(self.y_nominal, starImg.center.y, starImg.ymin)
            del psf
        return True, pos_shear

    def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
        disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0)
        disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
        disk = disk.shear(disk_shape)

        bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0)
        bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
        bulge = bulge.shear(bulge_shape)

        gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
        gal = gal.withFlux(flux)
        gal_shear = galsim.Shear(g1=g1, g2=g2)
        gal = gal.shear(gal_shear)
        final = galsim.Convolve(psf, gal)
        return final

    def rotateEllipticity(self, rotation):
        if rotation == 1:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e2_disk, self.e1_disk, -self.e2_bulge, self.e1_bulge, -self.e2_total, self.e1_total
        if rotation == 2:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e1_disk, -self.e2_disk, -self.e1_bulge, -self.e2_bulge, -self.e1_total, -self.e2_total
        if rotation == 3:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = self.e2_disk, -self.e1_disk, self.e2_bulge, -self.e1_bulge, self.e2_total, -self.e1_total

    def drawObject(self, img, final, noise_level=0.0, flux=None, filt=None, tel=None, exptime=150.):
        """ Override the method in parent class 
        Need to constrain the size of image stamp for extended objects
        """
        isUpdated = True
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
        stamp = final.drawImage(wcs=self.localWCS, offset=self.offset)
        stamp_arr = stamp.array
        mask = (stamp_arr >= 0.001*noise_level) # why 0.001?
        err = int(np.sqrt(mask.sum()))
        if np.mod(err, 2) == 1:
            err += 1
        # if err == 1:
        if err == 0:
            subSize = 16 # why 16?
        else:
            subSize = max([err, 16])
            fluxRatio = flux / stamp_arr[mask].sum()
            final = final.withScaledFlux(fluxRatio)

        imgSub = galsim.ImageF(subSize, subSize)

        # Draw with FFT
        # stamp = final.drawImage(image=imgSub, wcs=self.localWCS, offset=self.offset)

        # Draw with Photon Shoot
        stamp = final.drawImage(image=imgSub, wcs=self.localWCS, method='phot', offset=self.offset)
        
        stamp.setCenter(self.x_nominal, self.y_nominal)
        if np.sum(np.isnan(stamp.array)) >= 1:
            stamp.setZero()
        bounds = stamp.bounds & img.bounds
        if bounds.area() == 0:
            isUpdated = False
        else:
            img[bounds] += stamp[bounds]
        return img, stamp, isUpdated


    def getObservedEll(self, g1=0, g2=0):
        e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
        return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs