Chip.py 44.6 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
3
import galsim
import os
import numpy as np
Fang Yuedong's avatar
Fang Yuedong committed
4
5
import pickle
import json
Fang Yuedong's avatar
Fang Yuedong committed
6
7
from astropy.table import Table
from numpy.random import Generator, PCG64
Fang Yuedong's avatar
Fang Yuedong committed
8
9
from astropy.io import fits
from datetime import datetime
Fang Yuedong's avatar
Fang Yuedong committed
10

11
12
13
from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Instrument.FocalPlane import FocalPlane
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
Fang Yuedong's avatar
Fang Yuedong committed
14
from ObservationSim.Instrument._util import rotate_conterclockwise
15

Fang Yuedong's avatar
Fang Yuedong committed
16
17
18
19
20
21
try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources

Fang Yuedong's avatar
Fang Yuedong committed
22
class Chip(FocalPlane):
23
    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None):
Fang Yuedong's avatar
Fang Yuedong committed
24
25
26
        # Get focal plane (instance of paraent class) info
        # TODO: use chipID to config individual chip?
        super().__init__()
Fang Yuedong's avatar
Fang Yuedong committed
27
28
29
30
31
        # self.npix_x = 9216
        # self.npix_y = 9232
        # self.pix_scale  = 0.074 # pixel scale

        self._set_attributes_from_config(config)
Fang Yuedong's avatar
Fang Yuedong committed
32

33
34
        self.logger = logger

Fang Yuedong's avatar
Fang Yuedong committed
35
36
        # A chip ID must be assigned
        self.chipID = int(chipID)
Fang Yuedong's avatar
Fang Yuedong committed
37
        self.chip_name = str(chipID).rjust(2, '0')
Fang Yuedong's avatar
Fang Yuedong committed
38
39
40
41
42

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

Fang Yuedong's avatar
Fang Yuedong committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
        # [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":
60
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"]) == False:
Fang Yuedong's avatar
Fang Yuedong committed
61
62
63
64
65
66
67
68
69
70
71
72
73
                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
74
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"] == False):
Fang Yuedong's avatar
Fang Yuedong committed
75
76
77
78
79
80
81
82
83
84
85
                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)

Fang Yuedong's avatar
Fang Yuedong committed
86
87
88
89
90
91
        # 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:
92
93
94
95
96
97
            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)
Fang Yuedong's avatar
Fang Yuedong committed
98
        else:
Fang Yuedong's avatar
Fang Yuedong committed
99
100
            # self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
            self.sls_conf = []
101
102
103
104
105
106
107
108
109
110
111
112
            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))
Fang Yuedong's avatar
Fang Yuedong committed
113
114
115
116
        
        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

Fang Yuedong's avatar
Fang Yuedong committed
117
        # Define the sensor model
Fang Yuedong's avatar
Fang Yuedong committed
118
        if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
Fang Yuedong's avatar
Fang Yuedong committed
119
            self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
Fang Yuedong's avatar
Fang Yuedong committed
120
121
122
        else:
            self.sensor = galsim.Sensor()

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

Fang Yuedong's avatar
Fang Yuedong committed
125
126
127
128
129
130
131
132
133
134
135
    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])

Fang Yuedong's avatar
Fang Yuedong committed
136
137
138
139
140
141
142
143
144
145
146
    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"
Fang Yuedong's avatar
Fang Yuedong committed
147
        elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
Fang Yuedong's avatar
Fang Yuedong committed
148
            return "photometric"
Fang Yuedong's avatar
Fang Yuedong committed
149
150
        # elif self.filter_type in ["FGS"]:
        #     return "FGS"
Fang Yuedong's avatar
Fang Yuedong committed
151
152
153
154

    def _getChipEffCurve(self, filter_type):
        # CCD efficiency curves
        if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
Fang Yuedong's avatar
Fang Yuedong committed
155
        if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
Fang Yuedong's avatar
Fang Yuedong committed
156
157
        if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
        # Mirror efficiency:
158
159
160
161
        # 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
Fang Yuedong's avatar
Fang Yuedong committed
162
163
        # path = os.path.join(self.ccdEffCurve_dir, filename)
        # table = Table.read(path, format='ascii')
164
165
166
167
168
169
        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')
170
        throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
Fang Yuedong's avatar
Fang Yuedong committed
171
172
173
174
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
175
176
177
178
179
180
        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)
Fang Yuedong's avatar
Fang Yuedong committed
181
182
183
184

    def getChipFilter(self, chipID=None, filter_layout=None):
        """Return the filter index and type for a given chip #(chipID)
        """
Fang Yuedong's avatar
Fang Yuedong committed
185
        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
Fang Yuedong's avatar
Fang Yuedong committed
186
187
188
189
190
191
        if filter_layout is not None:
            return filter_layout[chipID][0], filter_layout[chipID][1]
        if chipID == None:
            chipID = self.chipID

        # updated configurations
Fang Yuedong's avatar
Fang Yuedong committed
192
193
        if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]")
        if chipID in [6, 15, 16, 25]:  filter_type = "y"
Fang Yuedong's avatar
Fang Yuedong committed
194
        if chipID in [11, 20]:         filter_type = "z"
Fang Yuedong's avatar
Fang Yuedong committed
195
        if chipID in [7, 24]:          filter_type = "i"
Fang Yuedong's avatar
Fang Yuedong committed
196
        if chipID in [14, 17]:         filter_type = "u"
Fang Yuedong's avatar
Fang Yuedong committed
197
198
199
200
201
202
203
        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'
Fang Yuedong's avatar
Fang Yuedong committed
204
205
206
207
208
209
        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
Fang Yuedong's avatar
Fang Yuedong committed
210
        NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
Fang Yuedong's avatar
Fang Yuedong committed
211
212
213
214
215
216
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
Fang Yuedong's avatar
Fang Yuedong committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
        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)
Fang Yuedong's avatar
Fang Yuedong committed
243
        else:
Fang Yuedong's avatar
Fang Yuedong committed
244
245
246
247
248
249
250
251
252
253
254
255
            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)
Fang Yuedong's avatar
Fang Yuedong committed
256
257
258
259



    def getSkyCoverage(self, wcs):
Fang Yuedong's avatar
Fang Yuedong committed
260
        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
Fang Yuedong's avatar
Fang Yuedong committed
261
262
263
264
265
266
267
268
269
270
        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)

Fang Yuedong's avatar
Fang Yuedong committed
271
    def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
Fang Yuedong's avatar
Fang Yuedong committed
272
        # magin in number of pix
Fang Yuedong's avatar
Fang Yuedong committed
273
274
275
276
277
278
279
280
        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")

Fang Yuedong's avatar
Fang Yuedong committed
281
282
        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
Fang Yuedong's avatar
Fang Yuedong committed
283
        if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
Fang Yuedong's avatar
Fang Yuedong committed
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
            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

Fang Yuedong's avatar
Fang Yuedong committed
307
308
309
310
311
312
313
    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, 
Fang Yuedong's avatar
Fang Yuedong committed
314
            pixel_scale=self.pix_scale,
Fang Yuedong's avatar
Fang Yuedong committed
315
316
317
            date=date_obs,
            time_obs=time_obs,
            im_type = im_type,
Fang Yuedong's avatar
Fang Yuedong committed
318
319
            exptime=exptime,
            chip_name=str(self.chipID).rjust(2, '0')
Fang Yuedong's avatar
Fang Yuedong committed
320
321
            )
        h_ext = generateExtensionHeader(
Fang Yuedong's avatar
Fang Yuedong committed
322
            chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
323
324
325
            xlen=self.npix_x, 
            ylen=self.npix_y, 
            ra=ra_cen, 
Fang Yuedong's avatar
Fang Yuedong committed
326
            dec=dec_cen,
Fang Yuedong's avatar
Fang Yuedong committed
327
328
329
330
331
            pa=img_rot.deg, 
            gain=self.gain, 
            readout=self.read_noise, 
            dark=self.dark_noise, 
            saturation=90000, 
Fang Yuedong's avatar
Fang Yuedong committed
332
333
334
335
            pixel_scale=self.pix_scale, 
            pixel_size=self.pix_size,
            xcen=self.x_cen,
            ycen=self.y_cen,
Fang Yuedong's avatar
Fang Yuedong committed
336
337
338
339
340
341
342
343
344
345
346
347
348
349
            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)
Fang Yuedong's avatar
Fang Yuedong committed
350
        hdu1.add_checksum()
Fang Yuedong's avatar
Fang Yuedong committed
351
        hdu2 = fits.ImageHDU(img.array, header=h_ext)
Fang Yuedong's avatar
Fang Yuedong committed
352
        hdu2.add_checksum()
Fang Yuedong's avatar
Fang Yuedong committed
353
354
355
356
        hdu1 = fits.HDUList([hdu1, hdu2])
        fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
        hdu1.writeto(fname, output_verify='ignore', overwrite=True)

357
    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):
Fang Yuedong's avatar
Fang Yuedong committed
358
359
360
361
362
363
        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"])
Fang Yuedong's avatar
Fang Yuedong committed
364
        fullwell = int(self.full_well)
Fang Yuedong's avatar
Fang Yuedong committed
365
        if config["ins_effects"]["add_hotpixels"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
366
367
368
            BoolHotPix = True
        else:
            BoolHotPix = False
Fang Yuedong's avatar
Fang Yuedong committed
369
        if config["ins_effects"]["add_deadpixels"]== True:
Fang Yuedong's avatar
Fang Yuedong committed
370
371
372
            BoolDeadPix = True
        else:
            BoolDeadPix = False
373
        self.logger = logger
Fang Yuedong's avatar
Fang Yuedong committed
374

375
376
377
378
379
        # 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.)

380
        # Add sky background
Zhang Xin's avatar
Zhang Xin committed
381
        if sky_map is None:
Fang Yuedong's avatar
Fang Yuedong committed
382
            sky_map = filt.getSkyNoise(exptime=exptime)
383
384
385
386
387
388
389
390
            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)
391
392
393
        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
Fang Yuedong's avatar
Fang Yuedong committed
394
            sky_map = sky_map * tel.pupil_area * exptime
Fang Yuedong's avatar
Fang Yuedong committed
395
        if config["ins_effects"]["add_back"] == True:
396
397
398
            img += sky_map
        del sky_map

Fang Yuedong's avatar
Fang Yuedong committed
399
        # Apply flat-field large scale structure for one chip
Fang Yuedong's avatar
Fang Yuedong committed
400
        if config["ins_effects"]["flat_fielding"] == True:
401
402
403
404
405
406
407
            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)
Fang Yuedong's avatar
Fang Yuedong committed
408
409
            flat_img = effects.MakeFlatSmooth(
                img.bounds, 
Fang Yuedong's avatar
Fang Yuedong committed
410
                int(config["random_seeds"]["seed_flat"]))
Fang Yuedong's avatar
Fang Yuedong committed
411
            flat_normal = flat_img / np.mean(flat_img.array)
412
413
            if self.survey_type == "photometric":
                img *= flat_normal
Fang Yuedong's avatar
Fang Yuedong committed
414
            del flat_normal
Fang Yuedong's avatar
Fang Yuedong committed
415
            if config["output_setting"]["flat_output"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
416
417
418
                del flat_img

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

433
434
435
436
437
        # 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)
Fang Yuedong's avatar
Fang Yuedong committed
438
439

        # Add cosmic-rays
Fang Yuedong's avatar
Fang Yuedong committed
440
        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS':
441
442
443
444
            if self.logger is not None:
                self.logger.info(("  Adding Cosmic-Ray"))
            else:
                print("  Adding Cosmic-Ray", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
445
            cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
446
                xLen=self.npix_x, yLen=self.npix_y, 
Fang Yuedong's avatar
Fang Yuedong committed
447
448
                exTime=exptime+0.5*self.readout_time, 
                cr_pixelRatio=0.003*(exptime+0.5*self.readout_time)/600.,
Fang Yuedong's avatar
Fang Yuedong committed
449
450
                gain=self.gain, 
                attachedSizes=self.attachedSizes,
Fang Yuedong's avatar
Fang Yuedong committed
451
                seed=SeedCosmicRay+pointing_ID*30+self.chipID)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
Fang Yuedong's avatar
Fang Yuedong committed
452
453
454
455
            img += cr_map
            cr_map[cr_map > 65535] = 65535
            cr_map[cr_map < 0] = 0
            crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
Fang Yuedong's avatar
Fang Yuedong committed
456
            del cr_map
Fang Yuedong's avatar
Fang Yuedong committed
457
            # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
Fang Yuedong's avatar
Fang Yuedong committed
458
            # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
Fang Yuedong's avatar
Fang Yuedong committed
459
            datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
460
461
462
463
464
465
466
467
468
469
470
471
            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,
Fang Yuedong's avatar
Fang Yuedong committed
472
                exptime=exptime)
Fang Yuedong's avatar
Fang Yuedong committed
473
474
            del crmap_gsimg

475
        # Apply PRNU effect and output PRNU flat file:
Fang Yuedong's avatar
Fang Yuedong committed
476
        if config["ins_effects"]["prnu_effect"] == True:
477
478
479
480
            if self.logger is not None:
                self.logger.info("  Applying PRNU effect")
            else:
                print("  Applying PRNU effect", flush=True)
481
482
483
484
            prnu_img = effects.PRNU_Img(
                xsize=self.npix_x, 
                ysize=self.npix_y, 
                sigma=0.01, 
Fang Yuedong's avatar
Fang Yuedong committed
485
                seed=int(config["random_seeds"]["seed_prnu"]+self.chipID))
486
            img *= prnu_img
Fang Yuedong's avatar
Fang Yuedong committed
487
            if config["output_setting"]["prnu_output"] == True:
488
                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
Fang Yuedong's avatar
Fang Yuedong committed
489
            if config["output_setting"]["flat_output"] == False:
490
491
492
                del prnu_img

        # Add dark current
Fang Yuedong's avatar
Fang Yuedong committed
493
        if config["ins_effects"]["add_dark"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
494
            dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
495
496
497
498
499
500
501
502
            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 
Fang Yuedong's avatar
Fang Yuedong committed
503
        if config["ins_effects"]["add_badcolumns"] == True:
504
            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
505

Fang Yuedong's avatar
Fang Yuedong committed
506
        # Add Bias level
Fang Yuedong's avatar
Fang Yuedong committed
507
        if config["ins_effects"]["add_bias"] == True:
508
509
510
511
            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")
512
513
            if config["ins_effects"]["bias_16channel"] == True:
                img = effects.AddBiasNonUniform16(img, 
Fang Yuedong's avatar
Fang Yuedong committed
514
                    bias_level=float(self.bias_level), 
515
516
517
518
519
                    nsecy = 2, nsecx=8, 
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
            elif config["ins_effects"]["bias_16channel"] == False:
                img += self.bias_level
Fang Yuedong's avatar
Fang Yuedong committed
520

521
        # Apply Nonlinearity on the chip image
Fang Yuedong's avatar
Fang Yuedong committed
522
        if config["ins_effects"]["non_linear"] == True:
523
524
525
526
            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)
527
528
529
            img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0)

        # Apply CCD Saturation & Blooming
Fang Yuedong's avatar
Fang Yuedong committed
530
        if config["ins_effects"]["saturbloom"] == True:
531
532
533
534
            if self.logger is not None:
                self.logger.info("  Applying CCD Saturation & Blooming")
            else:
                print("  Applying CCD Saturation & Blooming")
535
536
537
            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

        # Apply CTE Effect
Fang Yuedong's avatar
Fang Yuedong committed
538
        if config["ins_effects"]["cte_trail"] == True:
539
540
541
542
            if self.logger is not None:
                self.logger.info("  Apply CTE Effect")
            else:
                print("  Apply CTE Effect")
543
544
545
            img = effects.CTE_Effect(GSImage=img, threshold=27)

        # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
546
547
        if config["ins_effects"]["add_readout"] == True:
            seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID
548
549
550
551
552
            rng_readout = galsim.BaseDeviate(seed)
            readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
            img.addNoise(readout_noise)

        # Apply Gain & Quantization
553
554
555
556
        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)
557
558
559
560
561
562
563
564
565
        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
            
566
567
568
569
570
571
572
        img.array[img.array > 65535] = 65535
        img.replaceNegative(replace_value=0)
        img.quantize()

        ######################################################################################
        # Output images for calibration pointing
        ######################################################################################
Fang Yuedong's avatar
Fang Yuedong committed
573
        # Bias output
574
        if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
575
576
577
578
            if self.logger is not None:
                self.logger.info("  Output N frame Bias files")
            else:
                print("  Output N frame Bias files", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
579
            NBias = int(config["output_setting"]["NBias"])
Fang Yuedong's avatar
Fang Yuedong committed
580
581
582
            for i in range(NBias):
                BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                    self.npix_x, self.npix_y, 
Fang Yuedong's avatar
Fang Yuedong committed
583
                    bias_level=float(self.bias_level), 
Fang Yuedong's avatar
Fang Yuedong committed
584
                    ncombine=1, read_noise=self.read_noise, gain=1,
585
586
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
587
                # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced.
Fang Yuedong's avatar
Fang Yuedong committed
588
589
590
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
591
592
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=0.01, 
Fang Yuedong's avatar
Fang Yuedong committed
593
                            cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
594
595
596
597
598
599
600
                            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

Fang Yuedong's avatar
Fang Yuedong committed
601
                # Non-Linearity for Bias
Fang Yuedong's avatar
Fang Yuedong committed
602
                if config["ins_effects"]["non_linear"] == True:
603
604
605
606
                    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)
607
608
609
                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0)

                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
610
                if config["ins_effects"]["add_badcolumns"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
611
                    BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
Fang Yuedong's avatar
Fang Yuedong committed
612
613
614

                BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
615
616
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
617
618
                # BiasCombImg = effects.AddOverscan(
                #     BiasCombImg, 
Fang Yuedong's avatar
Fang Yuedong committed
619
                #     overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain, 
Fang Yuedong's avatar
Fang Yuedong committed
620
621
622
623
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                BiasCombImg.replaceNegative(replace_value=0)
                BiasCombImg.quantize()
                BiasCombImg = galsim.ImageUS(BiasCombImg)
Fang Yuedong's avatar
Fang Yuedong committed
624
                # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
Fang Yuedong's avatar
Fang Yuedong committed
625
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
626
627
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
628
                timestamp_obs += 10 * 60
Fang Yuedong's avatar
Fang Yuedong committed
629
630
631
632
633
                self.outputCal(
                    img=BiasCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
634
                    im_type='BIAS',
Fang Yuedong's avatar
Fang Yuedong committed
635
636
637
638
639
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=0.0)
Fang Yuedong's avatar
Fang Yuedong committed
640
641
642
            del BiasCombImg

        # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
643
        if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
644
645
646
647
            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)
Fang Yuedong's avatar
Fang Yuedong committed
648
            NFlat = int(config["output_setting"]["NFlat"])
Fang Yuedong's avatar
Fang Yuedong committed
649
            if config["ins_effects"]["add_bias"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
650
651
                biaslevel = self.bias_level
                overscan = biaslevel-2
652
            elif config["ins_effects"]["add_bias"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
653
654
                biaslevel = 0
                overscan = 0
655
            darklevel = self.dark_noise*(self.flat_exptime+0.5*self.readout_time)
Fang Yuedong's avatar
Fang Yuedong committed
656
657
658
659
660
661
662
663
            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, 
664
                    biaslevel=0,
665
666
                    seed_bias=SeedDefective+self.chipID,
                    logger=self.logger
Fang Yuedong's avatar
Fang Yuedong committed
667
                    )
Fang Yuedong's avatar
Fang Yuedong committed
668
669
670
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
671
                            xLen=self.npix_x, yLen=self.npix_y, 
672
                            exTime=self.flat_exptime+0.5*self.readout_time, 
Fang Yuedong's avatar
Fang Yuedong committed
673
                            cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
674
675
676
677
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+3)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
Fang Yuedong's avatar
Fang Yuedong committed
678
                    FlatCombImg += cr_map
Fang Yuedong's avatar
Fang Yuedong committed
679
                    del cr_map
Fang Yuedong's avatar
Fang Yuedong committed
680

Fang Yuedong's avatar
Fang Yuedong committed
681
                if config["ins_effects"]["non_linear"] == True:
682
683
684
685
                    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)
686
                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0)
Fang Yuedong's avatar
Fang Yuedong committed
687

Fang Yuedong's avatar
Fang Yuedong committed
688
                if config["ins_effects"]["cte_trail"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
689
690
                    FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)

691
692
693
694
695
                # 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)

Fang Yuedong's avatar
Fang Yuedong committed
696
                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
697
                if config["ins_effects"]["add_badcolumns"] == True:
698
                    FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
699

700
                # Add Bias level
Fang Yuedong's avatar
Fang Yuedong committed
701
                if config["ins_effects"]["add_bias"] == True:
702
703
704
705
                    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")
Fang Yuedong's avatar
Fang Yuedong committed
706
                    # img += float(config["ins_effects"]["bias_level"])
707
708
709
                    FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, 
                        bias_level=biaslevel, 
                        nsecy = 2, nsecx=8, 
710
711
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)
712
713
                
                # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
714
                if config["ins_effects"]["add_readout"] == True:
715
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 3
716
717
718
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    FlatCombImg.addNoise(readout_noise)
Fang Yuedong's avatar
Fang Yuedong committed
719
720
721

                FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
722
723
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
724
725
726
727
                # 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)
Fang Yuedong's avatar
Fang Yuedong committed
728
                # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
Fang Yuedong's avatar
Fang Yuedong committed
729
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
730
731
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
732
                timestamp_obs += 10 * 60
Fang Yuedong's avatar
Fang Yuedong committed
733
734
735
736
737
                self.outputCal(
                    img=FlatCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
738
                    im_type='FLAT',
Fang Yuedong's avatar
Fang Yuedong committed
739
740
741
742
743
744
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=self.flat_exptime)

Fang Yuedong's avatar
Fang Yuedong committed
745
746
747
748
749
750
751
            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
752
        if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
753
754
755
756
            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)
Fang Yuedong's avatar
Fang Yuedong committed
757
            NDark = int(config["output_setting"]["NDark"])
Fang Yuedong's avatar
Fang Yuedong committed
758
            if config["ins_effects"]["add_bias"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
759
760
                biaslevel = self.bias_level
                overscan = biaslevel-2
761
            elif config["ins_effects"]["add_bias"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
762
763
764
765
766
                biaslevel = 0
                overscan = 0
            for i in range(NDark):
                DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                    self.npix_x, self.npix_y, 
767
                    overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time,
Fang Yuedong's avatar
Fang Yuedong committed
768
                    ncombine=1, read_noise=self.read_noise, 
769
770
                    gain=1, seed_bias=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
771
772
773
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
774
                            xLen=self.npix_x, yLen=self.npix_y, 
775
                            exTime=self.dark_exptime+0.5*self.readout_time, 
Fang Yuedong's avatar
Fang Yuedong committed
776
                            cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
777
778
779
780
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+2)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
Fang Yuedong's avatar
Fang Yuedong committed
781
                    DarkCombImg += cr_map
Fang Yuedong's avatar
Fang Yuedong committed
782
783
784
785
                    cr_map[cr_map > 65535] = 65535
                    cr_map[cr_map < 0] = 0
                    crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
                    del cr_map
Fang Yuedong's avatar
Fang Yuedong committed
786
                    datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
787
788
789
790
791
792
793
794
795
796
797
798
799
800
                    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
Fang Yuedong's avatar
Fang Yuedong committed
801
802

                # Non-Linearity for Dark
Fang Yuedong's avatar
Fang Yuedong committed
803
                if config["ins_effects"]["non_linear"] == True:
804
805
806
807
                    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)
808
                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0)
Fang Yuedong's avatar
Fang Yuedong committed
809

Fang Yuedong's avatar
Fang Yuedong committed
810
                if config["ins_effects"]["cte_trail"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
811
812
                    DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)

813
814
815
816
817
                # 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)

Fang Yuedong's avatar
Fang Yuedong committed
818
                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
819
                if config["ins_effects"]["add_badcolumns"] == True:
820
                    DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
821

822
                # Add Bias level
Fang Yuedong's avatar
Fang Yuedong committed
823
                if config["ins_effects"]["add_bias"] == True:
824
825
826
827
                    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")
Fang Yuedong's avatar
Fang Yuedong committed
828
                    # img += float(config["ins_effects"]["bias_level"])
829
830
831
                    DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, 
                        bias_level=biaslevel, 
                        nsecy = 2, nsecx=8, 
832
833
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)
834
835

                # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
836
                if config["ins_effects"]["add_readout"] == True:
837
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 2
838
839
840
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    DarkCombImg.addNoise(readout_noise)
Fang Yuedong's avatar
Fang Yuedong committed
841
842
843
844

                DarkCombImg = effects.ApplyGainNonUniform16(
                    DarkCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
845
846
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
847
848
849
850
851
852
853
                # 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)
Fang Yuedong's avatar
Fang Yuedong committed
854
                # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
Fang Yuedong's avatar
Fang Yuedong committed
855
                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
856
857
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
858
                timestamp_obs += 10 * 60
Fang Yuedong's avatar
Fang Yuedong committed
859
860
861
862
863
                self.outputCal(
                    img=DarkCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
864
                    im_type='DARK',
Fang Yuedong's avatar
Fang Yuedong committed
865
866
867
868
869
                    pointing_ID=pointing_ID,
                    date_obs=date_obs,
                    time_obs=time_obs,
                    output_dir=chip_output.subdir,
                    exptime=self.dark_exptime)
Fang Yuedong's avatar
Fang Yuedong committed
870
871
872
873
            del DarkCombImg
        # img = galsim.ImageUS(img)

        # # 16 output channel, with each a single image file
Fang Yuedong's avatar
Fang Yuedong committed
874
        # if config["ins_effects"]["readout16"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
875
876
877
878
879
880
881
882
883
884
885
886
        #     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
Zhang Xin's avatar
Zhang Xin committed
887
        return img
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904

    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