Chip.py

import galsim
import os
import numpy as np
import pickle
import json
from astropy.table import Table
from numpy.random import Generator, PCG64
from astropy.io import fits
from datetime import datetime

from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Instrument.FocalPlane import FocalPlane
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument._util import rotate_conterclockwise

try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources

class Chip(FocalPlane):
    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None):
        # Get focal plane (instance of paraent class) info
        # TODO: use chipID to config individual chip?
        super().__init__()
        # self.npix_x = 9216
        # self.npix_y = 9232
        # self.pix_scale  = 0.074 # pixel scale

        self._set_attributes_from_config(config)

        self.logger = logger

        # A chip ID must be assigned
        self.chipID = int(chipID)
        self.chip_name = str(chipID).rjust(2, '0')

        # Get corresponding filter info
        self.filter_id, self.filter_type = self.getChipFilter()
        self.survey_type = self._getSurveyType()

        # [TODO]
        if self.filter_type != "FGS":
            self._getChipRowCol()

        # Set the relavent specs for FGS detectors
        # [TODO]
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        for key in chip_dict:
            setattr(self, key, chip_dict[key])
        if self.filter_type == "FGS":
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"]) == False:
                self.fdModel = None
            else:
                fgs_name = self.chip_name[0:4]
                try:
                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
                except AttributeError:
                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
        else:
            # Get the corresponding field distortion model
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"] == False):
                self.fdModel = None
            else:
                try:
                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
                except AttributeError:
                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)

        # Get boundary (in pix)
        self.bound = self.getChipLim()
        self.ccdEffCurve_dir = ccdEffCurve_dir
        self.CRdata_dir = CRdata_dir
        slsconfs = self.getChipSLSConf()
        if np.size(slsconfs) == 1:
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
        else:
            # self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
            self.sls_conf = []
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))
        
        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

        # Define the sensor model
        if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
            self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
        else:
            self.sensor = galsim.Sensor()

        self.flat_cube = None # for spectroscopic flat field cube simulation

    def _set_attributes_from_config(self, config):
        # Default setting
        self.read_noise = 5.0   # e/pix
        self.dark_noise = 0.02  # e/pix/s
        self.rotate_angle = 0.
        self.overscan   = 1000
        # Override default values
        for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]:
            if key in config["ins_effects"]:
                setattr(self, key, config["ins_effects"][key])

    def _getChipRowCol(self):
        self.rowID, self.colID = self.getChipRowCol(self.chipID)

    def getChipRowCol(self, chipID):
        rowID = ((chipID - 1) % 5) + 1
        colID = 6 - ((chipID - 1) // 5)
        return rowID, colID

    def _getSurveyType(self):
        if self.filter_type in ["GI", "GV", "GU"]:
            return "spectroscopic"
        elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
            return "photometric"
        # elif self.filter_type in ["FGS"]:
        #     return "FGS"

    def _getChipEffCurve(self, filter_type):
        # CCD efficiency curves
        if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
        if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
        if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
        # Mirror efficiency:
        # if filter_type == 'nuv': mirror_eff = 0.54
        # if filter_type == 'u': mirror_eff = 0.68
        # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
        # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
        # path = os.path.join(self.ccdEffCurve_dir, filename)
        # table = Table.read(path, format='ascii')
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
        throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
        try:
            with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)

    def getChipFilter(self, chipID=None, filter_layout=None):
        """Return the filter index and type for a given chip #(chipID)
        """
        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
        if filter_layout is not None:
            return filter_layout[chipID][0], filter_layout[chipID][1]
        if chipID == None:
            chipID = self.chipID

        # updated configurations
        if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]")
        if chipID in [6, 15, 16, 25]:  filter_type = "y"
        if chipID in [11, 20]:         filter_type = "z"
        if chipID in [7, 24]:          filter_type = "i"
        if chipID in [14, 17]:         filter_type = "u"
        if chipID in [9, 22]:          filter_type = "r"
        if chipID in [12, 13, 18, 19]: filter_type = "nuv"
        if chipID in [8, 23]:          filter_type = "g"
        if chipID in [1, 10, 21, 30]:  filter_type = "GI"
        if chipID in [2, 5, 26, 29]:   filter_type = "GV"
        if chipID in [3, 4, 27, 28]:   filter_type = "GU"
        if chipID in range(31, 43):    filter_type = 'FGS'
        filter_id = filter_type_list.index(filter_type)

        return filter_id, filter_type

    def getChipLim(self, chipID=None):
        """Calculate the edges in pixel for a given CCD chip on the focal plane
        NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
        if ((chipID is not None) and (int(chipID) <= 30)) or (self.chipID <= 30):
            # [TODO]
            if chipID == None:
                chipID = self.chipID
                rowID, colID = self.rowID, self.colID
            else:
                rowID, colID = self.getChipRowCol(chipID)
            gx1, gx2 = self.npix_gap_x
            gy = self.npix_gap_y

            # xlim of a given CCD chip
            xrem = 2*(colID - 1) - (self.nchip_x - 1)
            xcen = (self.npix_x//2 + gx1//2) * xrem
            if chipID >= 26 or chipID == 21:
                xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)
            if chipID <= 5 or chipID == 10:
                xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)
            nx0 = xcen - self.npix_x//2 + 1
            nx1 = xcen + self.npix_x//2

            # ylim of a given CCD chip
            yrem = (rowID - 1) - self.nchip_y // 2
            ycen = (self.npix_y + gy) * yrem
            ny0 = ycen - self.npix_y//2 + 1
            ny1 = ycen + self.npix_y//2
            return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1)
        else:
            xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10
            xcen = self.x_cen / self.pix_size
            ycen = self.y_cen / self.pix_size
            sign_x = [-1., 1., -1., 1.]
            sign_y = [-1., -1., 1., 1.]
            for i in range(4):
                x = xcen + sign_x[i] * self.npix_x / 2.
                y = ycen + sign_y[i] * self.npix_y / 2.
                x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
                xmin, xmax = min(xmin, x), max(xmax, x)
                ymin, ymax = min(ymin, y), max(ymax, y)
            return galsim.BoundsD(xmin, xmax, ymin, ymax)



    def getSkyCoverage(self, wcs):
        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
        return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)


    def getSkyCoverageEnlarged(self, wcs, margin=0.5):
        """The enlarged sky coverage of the chip
        """
        margin /= 60.0
        bound = self.getSkyCoverage(wcs)
        return galsim.BoundsD(bound.xmin - margin, bound.xmax + margin, bound.ymin - margin, bound.ymax + margin)

    def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
        # magin in number of pix
        if (ra_obj is not None) and (dec_obj is not None):
            if wcs is None:
                wcs = self.img.wcs
            pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
            x_image, y_image = pos_obj.x, pos_obj.y
        elif (x_image is None) or (y_image is None):
            raise ValueError("Either (ra_obj, dec_obj) or (x_image, y_image) should be given")

        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
        if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
            return False
        return True

    def getChipNoise(self, exptime=150.0):
        noise = self.dark_noise * exptime + self.read_noise**2
        return noise

    def getChipSLSConf(self):
        confFile = ''
        if self.chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
        if self.chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf']
        if self.chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf']
        if self.chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf']
        if self.chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf']
        if self.chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf']
        if self.chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf']
        if self.chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf']
        if self.chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf']
        if self.chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf']
        if self.chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf']
        if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
        return confFile

    def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, exptime=150.):
        h_prim = generatePrimaryHeader(
            xlen=self.npix_x, 
            ylen=self.npix_y, 
            pointNum = str(pointing_ID),
            ra=ra_cen, 
            dec=dec_cen, 
            pixel_scale=self.pix_scale,
            date=date_obs,
            time_obs=time_obs,
            im_type = im_type,
            exptime=exptime,
            chip_name=str(self.chipID).rjust(2, '0')
            )
        h_ext = generateExtensionHeader(
            chip=self,
            xlen=self.npix_x, 
            ylen=self.npix_y, 
            ra=ra_cen, 
            dec=dec_cen,
            pa=img_rot.deg, 
            gain=self.gain, 
            readout=self.read_noise, 
            dark=self.dark_noise, 
            saturation=90000, 
            pixel_scale=self.pix_scale, 
            pixel_size=self.pix_size,
            xcen=self.x_cen,
            ycen=self.y_cen,
            extName='raw')
        return h_prim, h_ext

    def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, output_dir, exptime=150.):
        h_prim, h_ext = self.generateHeader(
            ra_cen=ra_cen,
            dec_cen=dec_cen,
            img_rot=img_rot,
            im_type=im_type,
            pointing_ID=pointing_ID,
            date_obs=date_obs,
            time_obs=time_obs,
            exptime=exptime)
        hdu1 = fits.PrimaryHDU(header=h_prim)
        hdu1.add_checksum()
        hdu2 = fits.ImageHDU(img.array, header=h_ext)
        hdu2.add_checksum()
        hdu1 = fits.HDUList([hdu1, hdu2])
        fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
        hdu1.writeto(fname, output_verify='ignore', overwrite=True)

    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None):
        SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
        SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
        SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"])
        SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"])
        SeedDefective = int(config["random_seeds"]["seed_defective"])
        SeedCosmicRay = int(config["random_seeds"]["seed_CR"])
        fullwell = int(self.full_well)
        if config["ins_effects"]["add_hotpixels"] == True:
            BoolHotPix = True
        else:
            BoolHotPix = False
        if config["ins_effects"]["add_deadpixels"]== True:
            BoolDeadPix = True
        else:
            BoolDeadPix = False
        self.logger = logger

        # Get Poisson noise generator
        seed = int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID
        rng_poisson = galsim.BaseDeviate(seed)
        poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.)

        # Add sky background
        if sky_map is None:
            sky_map = filt.getSkyNoise(exptime=exptime)
            sky_map = sky_map * np.ones_like(img.array)
            sky_map = galsim.Image(array=sky_map)
            # Apply Poisson noise to the sky map
            # (NOTE): only for photometric chips
            # since it utilize the photon shooting
            # to draw stamps
            if self.survey_type == "photometric":
                sky_map.addNoise(poisson_noise)
        elif img.array.shape != sky_map.shape:
            raise ValueError("The shape img and sky_map must be equal.")
        elif tel is not None: # If sky_map is given in flux
            sky_map = sky_map * tel.pupil_area * exptime
        if config["ins_effects"]["add_back"] == True:
            img += sky_map
        del sky_map

        # Apply flat-field large scale structure for one chip
        if config["ins_effects"]["flat_fielding"] == True:
            if self.logger is not None:
                self.logger.info("  Creating and applying Flat-Fielding")
                msg = str(img.bounds)
                self.logger.info(msg)
            else:
                print("  Creating and applying Flat-Fielding", flush=True)
                print(img.bounds, flush=True)
            flat_img = effects.MakeFlatSmooth(
                img.bounds, 
                int(config["random_seeds"]["seed_flat"]))
            flat_normal = flat_img / np.mean(flat_img.array)
            if self.survey_type == "photometric":
                img *= flat_normal
            del flat_normal
            if config["output_setting"]["flat_output"] == False:
                del flat_img

        # Apply Shutter-effect for one chip
        if config["ins_effects"]["shutter_effect"] == True:
            if self.logger is not None:
                self.logger.info("  Apply shutter effect")
            else:
                print("  Apply shutter effect", flush=True)
            shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3)    # shutter effect normalized image for this chip
            if self.survey_type == "photometric":
                img *= shuttimg
            if config["output_setting"]["shutter_output"] == True:    # output 16-bit shutter effect image with pixel value <=65535
                shutt_gsimg = galsim.ImageUS(shuttimg*6E4)
                shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID))
                del shutt_gsimg
            del shuttimg

        # Add Poisson noise to the resulting images
        # (NOTE): this can only applied to the slitless image
        # since it dose not use photon shooting to draw stamps
        if self.survey_type == "spectroscopic":
            img.addNoise(poisson_noise)

        # Add cosmic-rays
        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS':
            if self.logger is not None:
                self.logger.info(("  Adding Cosmic-Ray"))
            else:
                print("  Adding Cosmic-Ray", flush=True)
            cr_map, cr_event_num = effects.produceCR_Map(
                xLen=self.npix_x, yLen=self.npix_y, 
                exTime=exptime+0.5*self.readout_time, 
                cr_pixelRatio=0.003*(exptime+0.5*self.readout_time)/600.,
                gain=self.gain, 
                attachedSizes=self.attachedSizes,
                seed=SeedCosmicRay+pointing_ID*30+self.chipID)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
            img += cr_map
            cr_map[cr_map > 65535] = 65535
            cr_map[cr_map < 0] = 0
            crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
            del cr_map
            # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
            # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
            datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
            date_obs = datetime_obs.strftime("%y%m%d")
            time_obs = datetime_obs.strftime("%H%M%S")
            self.outputCal(
                img=crmap_gsimg,
                ra_cen=ra_cen,
                dec_cen=dec_cen,
                img_rot=img_rot,
                im_type='CRS',
                pointing_ID=pointing_ID,
                date_obs=date_obs,
                time_obs=time_obs,
                output_dir=chip_output.subdir,
                exptime=exptime)
            del crmap_gsimg

        # Apply PRNU effect and output PRNU flat file:
        if config["ins_effects"]["prnu_effect"] == True:
            if self.logger is not None:
                self.logger.info("  Applying PRNU effect")
            else:
                print("  Applying PRNU effect", flush=True)
            prnu_img = effects.PRNU_Img(
                xsize=self.npix_x, 
                ysize=self.npix_y, 
                sigma=0.01, 
                seed=int(config["random_seeds"]["seed_prnu"]+self.chipID))
            img *= prnu_img
            if config["output_setting"]["prnu_output"] == True:
                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
            if config["output_setting"]["flat_output"] == False:
                del prnu_img

        # Add dark current
        if config["ins_effects"]["add_dark"] == True:
            dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
            img.addNoise(dark_noise)

        # Add Hot Pixels or/and Dead Pixels
        rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
        badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
        img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

        # Apply Bad lines 
        if config["ins_effects"]["add_badcolumns"] == True:
            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

        # Add Bias level
        if config["ins_effects"]["add_bias"] == True:
            if self.logger is not None:
                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
            else:
                print("  Adding Bias level and 16-channel non-uniformity")
            if config["ins_effects"]["bias_16channel"] == True:
                img = effects.AddBiasNonUniform16(img, 
                    bias_level=float(self.bias_level), 
                    nsecy = 2, nsecx=8, 
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
            elif config["ins_effects"]["bias_16channel"] == False:
                img += self.bias_level

        # Apply Nonlinearity on the chip image
        if config["ins_effects"]["non_linear"] == True:
            if self.logger is not None:
                self.logger.info("  Applying Non-Linearity on the chip image")
            else:
                print("  Applying Non-Linearity on the chip image", flush=True)
            img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0)

        # Apply CCD Saturation & Blooming
        if config["ins_effects"]["saturbloom"] == True:
            if self.logger is not None:
                self.logger.info("  Applying CCD Saturation & Blooming")
            else:
                print("  Applying CCD Saturation & Blooming")
            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

        # Apply CTE Effect
        if config["ins_effects"]["cte_trail"] == True:
            if self.logger is not None:
                self.logger.info("  Apply CTE Effect")
            else:
                print("  Apply CTE Effect")
            img = effects.CTE_Effect(GSImage=img, threshold=27)

        # Add Read-out Noise
        if config["ins_effects"]["add_readout"] == True:
            seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID
            rng_readout = galsim.BaseDeviate(seed)
            readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
            img.addNoise(readout_noise)

        # Apply Gain & Quantization
        if self.logger is not None:
            self.logger.info("  Applying Gain (and 16 channel non-uniformity) & Quantization")
        else:
            print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
        if config["ins_effects"]["gain_16channel"] == True:
            img = effects.ApplyGainNonUniform16(
                img, gain=self.gain, 
                nsecy = 2, nsecx=8, 
                seed=SeedGainNonuni+self.chipID,
                logger=self.logger)
        elif config["ins_effects"]["gain_16channel"] == False:
            img /= self.gain
            
        img.array[img.array > 65535] = 65535
        img.replaceNegative(replace_value=0)
        img.quantize()

        ######################################################################################
        # Output images for calibration pointing
        ######################################################################################
        # Bias output
        if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
            if self.logger is not None:
                self.logger.info("  Output N frame Bias files")
            else:
                print("  Output N frame Bias files", flush=True)
            NBias = int(config["output_setting"]["NBias"])
            for i in range(NBias):
                BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                    self.npix_x, self.npix_y, 
                    bias_level=float(self.bias_level), 
                    ncombine=1, read_noise=self.read_noise, gain=1,
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
                # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced.
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=0.01, 
                            cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/150., 
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+1)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
                    BiasCombImg += cr_map
                    del cr_map

                # Non-Linearity for Bias
                if config["ins_effects"]["non_linear"] == True:
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Bias image")
                    else:
                        print("  Applying Non-Linearity on the Bias image", flush=True)
                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0)

                # Apply Bad lines 
                if config["ins_effects"]["add_badcolumns"] == True:
                    BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5

                BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
                # BiasCombImg = effects.AddOverscan(
                #     BiasCombImg, 
                #     overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain, 
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                BiasCombImg.replaceNegative(replace_value=0)
                BiasCombImg.quantize()
                BiasCombImg = galsim.ImageUS(BiasCombImg)
                # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 10 * 60
                self.outputCal(
                    img=BiasCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='BIAS',
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=0.0)
            del BiasCombImg

        # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
        if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
            if self.logger is not None:
                self.logger.info("  Output N frame Flat-Field files")
            else:
                print("  Output N frame Flat-Field files", flush=True)
            NFlat = int(config["output_setting"]["NFlat"])
            if config["ins_effects"]["add_bias"] == True:
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["ins_effects"]["add_bias"] == False:
                biaslevel = 0
                overscan = 0
            darklevel = self.dark_noise*(self.flat_exptime+0.5*self.readout_time)
            for i in range(NFlat):
                FlatSingle = flat_img * prnu_img + darklevel
                FlatCombImg, FlatTag = effects.MakeFlatNcomb(
                    flat_single_image=FlatSingle, 
                    ncombine=1, 
                    read_noise=self.read_noise,
                    gain=1, 
                    overscan=overscan, 
                    biaslevel=0,
                    seed_bias=SeedDefective+self.chipID,
                    logger=self.logger
                    )
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=self.flat_exptime+0.5*self.readout_time, 
                            cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/150., 
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+3)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
                    FlatCombImg += cr_map
                    del cr_map

                if config["ins_effects"]["non_linear"] == True:
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Flat image")
                    else:
                        print("  Applying Non-Linearity on the Flat image", flush=True)
                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0)

                if config["ins_effects"]["cte_trail"] == True:
                    FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)

                # Add Hot Pixels or/and Dead Pixels
                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
                FlatCombImg = effects.DefectivePixels(FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

                # Apply Bad lines 
                if config["ins_effects"]["add_badcolumns"] == True:
                    FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

                # Add Bias level
                if config["ins_effects"]["add_bias"] == True:
                    if self.logger is not None:
                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
                    else:
                        print("  Adding Bias level and 16-channel non-uniformity")
                    # img += float(config["ins_effects"]["bias_level"])
                    FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, 
                        bias_level=biaslevel, 
                        nsecy = 2, nsecx=8, 
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)
                
                # Add Read-out Noise
                if config["ins_effects"]["add_readout"] == True:
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 3
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    FlatCombImg.addNoise(readout_noise)

                FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
                # FlatCombImg = effects.AddOverscan(FlatCombImg, overscan=overscan, gain=self.gain, widthl=27, widthr=27, widtht=8, widthb=8)
                FlatCombImg.replaceNegative(replace_value=0)
                FlatCombImg.quantize()
                FlatCombImg = galsim.ImageUS(FlatCombImg)
                # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 10 * 60
                self.outputCal(
                    img=FlatCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='FLAT',
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=self.flat_exptime)

            del FlatCombImg, FlatSingle, prnu_img
            # flat_img.replaceNegative(replace_value=0)
            # flat_img.quantize()
            # galsim.ImageUS(flat_img).write("%s/FlatImg_Vignette_%s.fits" % (chip_output.subdir, self.chipID))
            del flat_img

        # Export Dark current images
        if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
            if self.logger is not None:
                self.logger.info("  Output N frame Dark Current files")
            else:
                print("  Output N frame Dark Current files", flush=True)
            NDark = int(config["output_setting"]["NDark"])
            if config["ins_effects"]["add_bias"] == True:
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["ins_effects"]["add_bias"] == False:
                biaslevel = 0
                overscan = 0
            for i in range(NDark):
                DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                    self.npix_x, self.npix_y, 
                    overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time,
                    ncombine=1, read_noise=self.read_noise, 
                    gain=1, seed_bias=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=self.dark_exptime+0.5*self.readout_time, 
                            cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/150., 
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+2)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
                    DarkCombImg += cr_map
                    cr_map[cr_map > 65535] = 65535
                    cr_map[cr_map < 0] = 0
                    crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
                    del cr_map
                    datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                    date_obs = datetime_obs.strftime("%y%m%d")
                    time_obs = datetime_obs.strftime("%H%M%S")
                    self.outputCal(
                        img=crmap_gsimg,
                        ra_cen=ra_cen,
                        dec_cen=dec_cen,
                        img_rot=img_rot,
                        im_type='CRD',
                        pointing_ID=pointing_ID,
                        date_obs=date_obs,
                        time_obs=time_obs,
                        output_dir=chip_output.subdir,
                        exptime=self.dark_exptime)
                    del crmap_gsimg

                # Non-Linearity for Dark
                if config["ins_effects"]["non_linear"] == True:
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Dark image")
                    else:
                        print("  Applying Non-Linearity on the Dark image", flush=True)
                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0)

                if config["ins_effects"]["cte_trail"] == True:
                    DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)

                # Add Hot Pixels or/and Dead Pixels
                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
                DarkCombImg = effects.DefectivePixels(DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

                # Apply Bad lines 
                if config["ins_effects"]["add_badcolumns"] == True:
                    DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

                # Add Bias level
                if config["ins_effects"]["add_bias"] == True:
                    if self.logger is not None:
                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
                    else:
                        print("  Adding Bias level and 16-channel non-uniformity")
                    # img += float(config["ins_effects"]["bias_level"])
                    DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, 
                        bias_level=biaslevel, 
                        nsecy = 2, nsecx=8, 
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)

                # Add Read-out Noise
                if config["ins_effects"]["add_readout"] == True:
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 2
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    DarkCombImg.addNoise(readout_noise)

                DarkCombImg = effects.ApplyGainNonUniform16(
                    DarkCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
                # DarkCombImg = effects.AddOverscan(
                #     DarkCombImg, 
                #     overscan=overscan, gain=self.gain, 
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                DarkCombImg.replaceNegative(replace_value=0)
                DarkCombImg.quantize()
                DarkCombImg = galsim.ImageUS(DarkCombImg)
                # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 10 * 60
                self.outputCal(
                    img=DarkCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='DARK',
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=self.dark_exptime)
            del DarkCombImg
        # img = galsim.ImageUS(img)

        # # 16 output channel, with each a single image file
        # if config["ins_effects"]["readout16"] == True:
        #     print("  16 Output Channel simulation")
        #     for coli in [0, 1]:
        #         for rowi in range(8):
        #             sub_img = effects.readout16(
        #                 GSImage=img, 
        #                 rowi=rowi, 
        #                 coli=coli, 
        #                 overscan_value=self.overscan)
        #             rowcoltag = str(rowi) + str(coli)
        #             img_name_root = chip_output.img_name.split(".")[0]
        #             sub_img.write("%s/%s_%s.fits" % (chip_output.subdir, img_name_root, rowcoltag))
        #     del sub_img
        return img

    def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'):
        from astropy.io import fits
        try:
            with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path:
                flat_fits = fits.open(data_path, ignore_missing_simple=True)
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path:
                flat_fits = fits.open(data_path, ignore_missing_simple=True)

        fl = len(flat_fits)
        fl_sh = flat_fits[0].data.shape
        assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!'
        self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]])
        for i in np.arange(0, fl, 1):
            self.flat_cube[i] = flat_fits[i].data