import os import galsim import numpy as np from astropy.io import fits from observation_sim.instruments.chip import chip_utils from observation_sim.instruments.chip import effects from astropy.time import Time from datetime import datetime, timezone def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param): self.chip_output.Log_info("Apply pre/over-scan") chip.img = chip_utils.AddPreScan(GSImage=chip.img, pre1=chip.prescan_x, pre2=chip.prescan_y, over1=chip.overscan_x, over2=chip.overscan_y) if obs_param["add_dark"] is True: ny = int(chip.npix_y/2) base_dark = (ny-1)*(chip.readout_time/ny)*chip.dark_noise chip.img.array[(chip.prescan_y+ny):-(chip.prescan_y+ny), :] = base_dark return chip, filt, tel, pointing def add_crosstalk(self, chip, filt, tel, pointing, catalog, obs_param): crosstalk=np.zeros([16, 16]) crosstalk[0, :] = np.array([1., 1e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[1, :] = np.array([1e-4, 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[2, :] = np.array([0., 0., 1., 1e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[3, :] = np.array([0., 0., 1e-4, 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[4, :] = np.array([0., 0., 0., 0., 1., 1e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[5, :] = np.array([0., 0., 0., 0., 1e-4, 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[6, :] = np.array([0., 0., 0., 0., 0., 0., 1., 1e-4, 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[7, :] = np.array([0., 0., 0., 0., 0., 0., 1e-4, 1., 0., 0., 0., 0., 0., 0., 0., 0.]) crosstalk[8, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 1., 1e-4, 0., 0., 0., 0., 0., 0.]) crosstalk[9, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 1e-4, 1., 0., 0., 0., 0., 0., 0.]) crosstalk[10, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1e-4, 0., 0., 0., 0.]) crosstalk[11, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1e-4, 1., 0., 0., 0., 0.]) crosstalk[12, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1e-4, 0., 0.]) crosstalk[13, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1e-4, 1., 0., 0.]) crosstalk[14, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1e-4]) crosstalk[15, :] = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1e-4, 1.]) # 2*8 -> 1*16 img = formatOutput(chip.img) ny, nx = img.array.shape nsecy = 1 nsecx = 16 dy = int(ny/nsecy) dx = int(nx/nsecx) newimg = galsim.Image(nx, ny, init_value=0) for i in range(16): for j in range(16): newimg.array[:, int(i*dx):int(i*dx+dx)] += crosstalk[i, j]*img.array[:, int(j*dx):int(j*dx+dx)] # 1*16 -> 2*8 newimg = formatRevert(newimg) chip.img.array[:, :] = newimg.array[:, :] return chip, filt, tel, pointing def add_readout_noise(self, chip, filt, tel, pointing, catalog, obs_param): seed = int(self.overall_config["random_seeds"] ["seed_readout"]) + pointing.id*30 + chip.chipID rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise( rng=rng_readout, sigma=chip.read_noise) chip.img.addNoise(readout_noise) return chip, filt, tel, pointing def apply_gain(self, chip, filt, tel, pointing, catalog, obs_param): self.chip_output.Log_info(" Applying Gain") if obs_param["gain_16channel"] is True: chip.img, chip.gain_channel = effects.ApplyGainNonUniform16(chip.img, gain=chip.gain, nsecy=chip.nsecy, nsecx=chip.nsecx, seed=self.overall_config["random_seeds"]["seed_gainNonUniform"]+chip.chipID) elif obs_param["gain_16channel"] is False: chip.img /= chip.gain chip.gain_channel = np.ones(chip.nsecy*chip.nsecx)*chip.gain return chip, filt, tel, pointing def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param): if not hasattr(self, 'h_ext'): _, _ = self.prepare_headers(chip=chip, pointing=pointing) self.updateHeaderInfo(header_flag='ext', keys=['SHTSTAT', 'SHTOPEN1', 'SHTCLOS0', 'SHTCLOS1', 'EXPTIME'], values=[ False, self.h_ext['SHTOPEN0'], self.h_ext['SHTOPEN0'], self.h_ext['SHTOPEN0'], 0.0]) # renew header info datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) t_obs = Time(datetime_obs) # ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光 t_obs_renew = Time(t_obs.mjd - 2. / 86400., format="mjd") t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp( t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) self.updateHeaderInfo(header_flag='prim', keys=[ 'DATE-OBS'], values=[t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) gains1 = list(chip.gain_channel[0:8]) gains2 = list(chip.gain_channel[8:]) gains2.reverse() gains = np.append(gains1, gains2) self.updateHeaderInfo(header_flag='ext', keys=['GAIN01', 'GAIN02', 'GAIN03', 'GAIN04', 'GAIN05', 'GAIN06', 'GAIN07', 'GAIN08', 'GAIN09', 'GAIN10', 'GAIN11', 'GAIN12', 'GAIN13', 'GAIN14', 'GAIN15', 'GAIN16'], values=gains) if obs_param["format_output"] is True: self.chip_output.Log_info(" Apply 1*16 format") chip.img = chip_utils.formatOutput(GSImage=chip.img) chip.nsecy = 1 chip.nsecx = 16 chip.img.array[chip.img.array > 65535] = 65535 chip.img.replaceNegative(replace_value=0) chip.img.quantize() chip.img = galsim.Image(chip.img.array, dtype=np.uint16) fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits') # f_name_size = 68 # if (len(self.h_prim['FILENAME']) > f_name_size): # self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[ # self.h_prim['FILENAME'][0:f_name_size]]) hdu1 = fits.PrimaryHDU(header=self.h_prim) self.updateHeaderInfo(header_flag='ext', keys=['DATASECT'], values=[ str(chip.img.array.shape[1]) + 'x' + str(chip.img.array.shape[0])]) hdu2 = fits.ImageHDU(chip.img.array, header=self.h_ext) hdu2.header.comments["XTENSION"] = "image extension" hdu = fits.HDUList([hdu1, hdu2]) hdu[0].add_datasum(when='data unit checksum') hdu[0].add_checksum(when='HDU checksum', override_datasum=True) hdu[1].add_datasum(when='data unit checksum') hdu[1].add_checksum(when='HDU checksum', override_datasum=True) hdu.writeto(fname, output_verify='ignore', overwrite=True) return chip, filt, tel, pointing