Chip.py

import galsim
import os
import numpy as np
import Instrument.Chip.Effects as effects
from Instrument.FocalPlane import FocalPlane
from astropy.table import Table
from numpy.random import Generator, PCG64
from Config.Header import generatePrimaryHeader, generateExtensionHeader
from astropy.io import fits
from datetime import datetime

class Chip(FocalPlane):
    def __init__(self, chipID, ccdEffCurve_dir, CRdata_dir, normalize_dir=None, sls_dir=None, config=None, treering_func=None):
        # Get focal plane (instance of paraent class) info
        # TODO: use chipID to config individual chip?
        super().__init__()
        # if config is not None:
        #     self.npix_x     = config["npix_x"]
        #     self.npix_y     = config["npix_y"]
        #     self.read_noise = config["read_noise"]
        #     self.dark_noise = config["dark_noise"]
        #     self.pix_scale  = config["pix_scale"]
        #     self.gain       = config["gain"]
        #     self.bias_level = config["bias_level"]
        #     self.overscan   = config["overscan"]
        # else:
        # Default setting
        self.npix_x = 9216
        self.npix_y = 9232
        self.read_noise = 5.0   # e/pix
        self.dark_noise = 0.02  # e/pix/s
        self.pix_scale  = 0.074 # pixel scale
        # self.gain       = 1.0
        # self.bias_level = 1000  # e-/pix
        self.gain = float(config["gain"])
        self.bias_level = float(config["bias_level"])
        self.overscan   = 1000
        self.exptime    = 150   # second
        # self.dark_exptime = 300  # second
        # self.flat_exptime = 150  # second
        self.dark_exptime = float(config['dark_exptime'])
        self.flat_exptime = float(config['flat_exptime'])

        # A chip ID must be assigned
        self.chipID = int(chipID)
        self._getChipRowCol()

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

        # Get boundary (in pix)
        self.bound = self.getChipLim()
        self.ccdEffCurve_dir = ccdEffCurve_dir
        self.CRdata_dir = CRdata_dir
        self.normalize_dir = normalize_dir
        self.sls_dir=sls_dir
        # self.sls_conf = os.path.join(self.sls_dir, self.getChipSLSConf())
        slsconfs = self.getChipSLSConf()
        if np.size(slsconfs) == 1:
            self.sls_conf = [os.path.join(self.sls_dir, slsconfs)]
        else:
            self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
        
        if self.normalize_dir is not None:
            self._getNormF()
        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

        # Define the sensor
        if config["bright_fatter"].lower() == "y":
            self.sensor = galsim.SiliconSensor(strength=config["df_strength"], treering_func=treering_func)
        else:
            self.sensor = galsim.Sensor()

    # def _getChipRowCol(self):
    #     self.rowID = (self.chipID - 1) // 5 + 1
    #     self.colID = (self.chipID - 1) % 5 + 1
    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"
        else:
            return "photometric"

    def _getNormF(self):
        self.normF_star = Table.read(os.path.join(self.normalize_dir, 'SLOAN_SDSS.g.fits'))
        self.normF_galaxy = Table.read(os.path.join(self.normalize_dir, 'lsst_throuput_g.fits'))

    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']: filename = 'Astro_MB.txt'
        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')
        throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
        path = os.path.join(self.CRdata_dir, 'wfc-cr-attachpixel.dat')
        self.attachedSizes = np.loadtxt(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"]
        # TODO: maybe a more elegent way other than hard coded?
        # e.g. use something like a nested dict:
        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>30 or chipID<1: raise ValueError("!!! Chip ID: [1,30]")
        # if chipID in [10, 15, 16, 21]: filter_type = 'y'
        # if chipID in [11, 20]:         filter_type = "z"
        # if chipID in [9, 22]:           filter_type = "i"
        # if chipID in [12, 19]:         filter_type = "u"
        # if chipID in [7, 24]:         filter_type = "r"
        # if chipID in [14, 13, 18, 17]:    filter_type = "nuv"
        # if chipID in [8, 23]:         filter_type = "g"
        # if chipID in [6, 5, 25, 26]:    filter_type = "GI"
        # if chipID in [27, 30, 1, 4]:    filter_type = "GV"
        # if chipID in [28, 29, 2, 3]:    filter_type = "GU"
        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"
        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 observation.
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
        # if chipID == None:
        #     chipID = self.chipID
        
        # gx = self.npix_gap_x
        # gy1, gy2 = self.npix_gap_y

        # # xlim of a given ccd chip
        # xrem = (chipID-1)%self.nchip_x - self.nchip_x // 2
        # xcen = (self.npix_x + gx) * xrem
        # nx0 = xcen - self.npix_x//2 + 1
        # nx1 = xcen + self.npix_x//2

        # # ylim of a given ccd chip
        # yrem = 2*((chipID-1)//self.nchip_x) - (self.nchip_y-1)
        # ycen = (self.npix_y//2 + gy1//2) * yrem
        # if chipID <= 6: ycen = (self.npix_y//2 + gy1//2) * yrem - (gy2-gy1)
        # if chipID >= 25: ycen = (self.npix_y//2 + gy1//2) * yrem + (gy2-gy1)
        # ny0 = ycen - self.npix_y//2 + 1
        # ny1 = ycen + self.npix_y//2

        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)


    def getSkyCoverage(self, wcs):
        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, dec_obj, wcs=None, margin=1):
        # magin in number of pix
        if wcs is None:
            wcs = self.img.wcs
        pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
        if (pos_obj.x - xmin) * (pos_obj.x - xmax) > 0.0 or (pos_obj.y - ymin) * (pos_obj.y - 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 addNoise_phot(self, img, exptime=150.0, sky_noise=0., seed=31415):
        rng = galsim.BaseDeviate(seed)
        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
        img.addNoise(dark_noise)
        ccd_noise = galsim.CCDNoise(rng, sky_level=sky_noise, gain=self.gain, read_noise=self.read_noise)
        img.addNoise(ccd_noise)
        return img

    def addNoise_spec(self, config, tel, img, sky_map, exptime=150.0, seed=31415):
        if img.array.shape != sky_map.shape:
            raise ValueError("The shape img and sky_map must be equal.")

        # n_img_arrar = (img.array + sky_map) * tel.pupil_area * exptime
        # Should be the following?
        rng = galsim.BaseDeviate(seed)
        noise_img = galsim.Image(img, copy=True)
        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
        img.addNoise(dark_noise)
        if config["abs_back"].lower() == "y":
            img += (sky_map * tel.pupil_area * exptime)
        ccd_noise = galsim.CCDNoise(rng, gain=self.gain, read_noise=self.read_noise)
        img.addNoise(ccd_noise)
        return img
        # pNoise = self.dark_noise * exptime
        # noise_img = galsim.Image(n_img_arrar, copy=True)
        # dNose = galsim.PoissonNoise(rng=rng, sky_level=pNoise)
        # noise_img.addNoise(dNose)
        # rdNoise = galsim.GaussianNoise(rng=rng, sigma=self.read_noise)
        # noise_img.addNoise(rdNoise)
        # return noise_img

    def addNoise(self, config, tel, filt, img, sky_map=None, exptime=150.0):
        if self.survey_type == "photometric":
            sky_level = filt.getSkyNoise(exptime=exptime, gain=self.gain)
            img = self.addNoise_phot(
                img=img, 
                exptime=exptime, 
                sky_noise=filt.getSkyNoise(exptime=exptime, gain=self.gain), 
                seed=int(config["seed_skynoise"]))
            if config["abs_back"].lower() == "y":
                img += sky_level
        elif self.survey_type == "spectroscopic":
            if sky_map is None:
                raise ValueError("Must provide sky_map for spectroscopic chip")
            img = self.addNoise_spec(
                config=config, 
                tel=tel, 
                img=img, 
                sky_map=sky_map, 
                exptime=exptime, 
                seed=int(config["seed_skynoise"]))
        return img

    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, 
            psize=self.pix_scale, 
            row_num=self.rowID, 
            col_num=self.colID,
            date=date_obs,
            time_obs=time_obs,
            im_type = im_type,
            exptime=exptime
            )
        h_ext = generateExtensionHeader(
            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, 
            psize=self.pix_scale, 
            row_num=self.rowID, 
            col_num=self.colID,
            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)
        hdu2 = fits.ImageHDU(img.array, header=h_ext)
        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'):
        SeedGainNonuni=int(config["seed_gainNonUniform"])
        SeedBiasNonuni=int(config["seed_biasNonUniform"])
        SeedRnNonuni = int(config["seed_rnNonUniform"])
        SeedBadColumns = int(config["seed_badcolumns"])
        SeedDefective = int(config["seed_defective"])
        SeedCosmicRay = int(config["seed_CR"])
        fullwell = int(config["full_well"])
        if config["add_hotpixels"].lower() == "y":
            BoolHotPix = True
        else:
            BoolHotPix = False
        if config["add_deadpixels"].lower() == "y":
            BoolDeadPix = True
        else:
            BoolDeadPix = False

        # Apply flat-field large scale structure for one chip
        if config["flat_fielding"].lower() == "y":
            print("  Creating and applying Flat-Fielding", flush=True)
            print(img.bounds, flush=True)
            flat_img = effects.MakeFlatSmooth(
                img.bounds, 
                int(config["seed_flat"]))
            flat_normal = flat_img / np.mean(flat_img.array)
            img *= flat_normal
            del flat_normal
            if config["flat_output"].lower() == "n":
                del flat_img

        # Apply Shutter-effect for one chip
        if config["shutter_effect"].lower() == "y":
            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
            img *= shuttimg
            if config["shutter_output"].lower() == "y":    # 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

        # Apply PRNU effect and output PRNU flat file:
        if config["prnu_effect"].lower() == "y":
            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["seed_prnu"]+self.chipID))
            img *= prnu_img
            if config["prnu_output"].lower() == "y":
                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
            if config["flat_output"].lower() == "n":
                del prnu_img

        # Add dark current
        img += self.dark_noise*exptime

        # Add cosmic-rays
        if config["cosmic_ray"].lower() == "y" and pointing_type=='MS':
            print("  Adding Cosmic-Ray", flush=True)
            cr_map = effects.produceCR_Map(
                xLen=self.npix_x, yLen=self.npix_y, 
                exTime=exptime, 
                cr_pixelRatio=0.003, 
                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.fromtimestamp(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=150.)
            del crmap_gsimg

        # Add Bias level
        if config["add_bias"].lower() == "y":
            print("  Adding Bias level and 16-channel non-uniformity")
            # img += float(config["bias_level"])
            img = effects.AddBiasNonUniform16(img, 
                bias_level=float(config["bias_level"]), 
                nsecy = 2, nsecx=8, 
                seed=SeedBiasNonuni+self.chipID)

        # Bias output
        if config["bias_output"].lower() == "y" and pointing_type=='CAL':
            print("  Output N frame Bias files", flush=True)
            NBias = int(config["NBias"])
            for i in range(NBias):
                BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                    self.npix_x, self.npix_y, 
                    bias_level=float(config["bias_level"]), 
                    ncombine=1, read_noise=self.read_noise, gain=1,
                    seed=SeedBiasNonuni+self.chipID)
                if config["cosmic_ray"].lower() == "y":
                    if config["cray_differ"].lower() == "y":
                        cr_map = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=0.01, 
                            cr_pixelRatio=0.003*0.01/self.exptime, 
                            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["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Bias image", flush=True)
                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=1.e-7, beta2=1.e-10)

                BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # BiasCombImg = effects.AddOverscan(
                #     BiasCombImg, 
                #     overscan=float(config["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.fromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 5 * 60
                self.outputCal(
                    img=BiasCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='CLB',
                    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["flat_output"].lower() == "y" and pointing_type=='CAL':
            print("  Output N frame Flat-Field files", flush=True)
            NFlat = int(config["NFlat"])
            if config["add_bias"].lower() == "y":
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["add_bias"].lower() == "n":
                biaslevel = 0
                overscan = 0
            darklevel = self.dark_noise*self.flat_exptime
            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=biaslevel,
                    seed_bias=SeedDefective+self.chipID
                    )
                if config["cosmic_ray"].lower() == "y":
                    if config["cray_differ"].lower() == "y":
                        cr_map = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=self.flat_exptime, 
                            cr_pixelRatio=0.003*self.flat_exptime/self.exptime, 
                            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["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Flat image", flush=True)
                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=1.e-7, beta2=1.e-10)

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

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

                # 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=self.bias_level)

                FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # 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.fromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 5 * 60
                self.outputCal(
                    img=FlatCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='CLF',
                    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["dark_output"].lower() == "y" and pointing_type=='CAL':
            print("  Output N frame Dark Current files", flush=True)
            NDark = int(config["NDark"])
            if config["add_bias"].lower() == "y":
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["add_bias"].lower() == "n":
                biaslevel = 0
                overscan = 0
            for i in range(NDark):
                DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                    self.npix_x, self.npix_y, 
                    overscan=overscan, bias_level=biaslevel, darkpsec=0.02, exptime=self.dark_exptime,
                    ncombine=1, read_noise=self.read_noise, 
                    gain=1, seed_bias=SeedBiasNonuni+self.chipID)
                if config["cosmic_ray"].lower() == "y":
                    if config["cray_differ"].lower() == "y":
                        cr_map = effects.produceCR_Map(
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=self.dark_exptime, 
                            cr_pixelRatio=0.003*self.dark_exptime/self.exptime, 
                            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.fromtimestamp(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["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Dark image", flush=True)
                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=1.e-7, beta2=1.e-10)

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

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

                # 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=self.bias_level)

                DarkCombImg = effects.ApplyGainNonUniform16(
                    DarkCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # 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.fromtimestamp(timestamp_obs)
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                timestamp_obs += 5 * 60
                self.outputCal(
                    img=DarkCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
                    im_type='CLD',
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=self.dark_exptime)
            del DarkCombImg

        # Apply Nonlinearity on the chip image
        if config["non_linear"].lower() == "y":
            print("  Applying Non-Linearity on the chip image", flush=True)
            img = effects.NonLinearity(GSImage=img, beta1=1.e-7, beta2=1.e-10)

        # Apply CCD Saturation & Blooming
        if config["saturbloom"].lower() == "y":
            print("  Applying CCD Saturation & Blooming")
            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

        # Apply CTE Effect
        if config["cte_trail"].lower() == "y":
            print("  Apply CTE Effect")
            img = effects.CTE_Effect(GSImage=img, threshold=27)

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

        # 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=self.bias_level)
        

        # Apply Gain & Quantization
        print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
        img = effects.ApplyGainNonUniform16(
            img, gain=self.gain, 
            nsecy = 2, nsecx=8, 
            seed=SeedGainNonuni+self.chipID)
        img.array[img.array > 65535] = 65535
        img.replaceNegative(replace_value=0)
        img.quantize()
        # img = galsim.ImageUS(img)

        # # 16 output channel, with each a single image file
        # if config["readout16"].lower() == "y":
        #     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