From a6664237cca40da77e9e100201f8fe8a1611037a Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Mon, 14 Aug 2023 05:58:47 +0800 Subject: [PATCH] new sim: modified ways to add detector effects: e.g. brighter-fatter etc. --- ObservationSim/Instrument/Chip/Chip.py | 345 ++------ ObservationSim/Instrument/Chip/ChipUtils.py | 196 +++++ .../Instrument/Chip/lib_bf/diffusion_X1.c | 120 +++ .../Instrument/Chip/lib_bf/nrutil.c | 769 ++++++++++++++++++ .../Instrument/Chip/lib_bf/nrutil.h | 94 +++ ObservationSim/MockObject/Galaxy.py | 130 +-- ObservationSim/MockObject/MockObject.py | 133 +-- ObservationSim/ObservationSim.py | 4 +- config/config_50sqdeg.yaml | 10 +- libmoduleDfBF/diffusion_X1.c | 120 +++ profile_C6.sh | 13 +- run_C6.pbs | 4 +- setup.py | 30 +- 13 files changed, 1491 insertions(+), 477 deletions(-) create mode 100644 ObservationSim/Instrument/Chip/ChipUtils.py create mode 100644 ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c create mode 100644 ObservationSim/Instrument/Chip/lib_bf/nrutil.c create mode 100644 ObservationSim/Instrument/Chip/lib_bf/nrutil.h create mode 100644 libmoduleDfBF/diffusion_X1.c diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index cc0b4ce..bb93b0a 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -12,6 +12,7 @@ from ObservationSim.Instrument.Chip import Effects as effects from ObservationSim.Instrument.FocalPlane import FocalPlane from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader from ObservationSim.Instrument._util import rotate_conterclockwise +from ObservationSim.Instrument.Chip import ChipUtils as chip_utils try: import importlib.resources as pkg_resources @@ -22,11 +23,7 @@ except ImportError: class Chip(FocalPlane): def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None): # Get focal plane (instance of paraent class) info - # TODO: use chipID to config individual chip? super().__init__() - # self.npix_x = 9216 - # self.npix_y = 9232 - # self.pix_scale = 0.074 # pixel scale self.nsecy = 2 self.nsecx = 8 self.gain_channel = np.ones(self.nsecy* self.nsecx) @@ -42,12 +39,10 @@ class Chip(FocalPlane): self.filter_id, self.filter_type = self.getChipFilter() self.survey_type = self._getSurveyType() - # [TODO] if self.filter_type != "FGS": self._getChipRowCol() - # Set the relavent specs for FGS detectors - # [TODO] + # Set the relavent specs for detectors try: with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition: with open(chip_definition, "r") as f: @@ -58,6 +53,7 @@ class Chip(FocalPlane): chip_dict = json.load(f)[str(self.chipID)] for key in chip_dict: setattr(self, key, chip_dict[key]) + if self.filter_type == "FGS": if ("field_dist" in config) and (config["ins_effects"]["field_dist"]) == False: self.fdModel = None @@ -77,7 +73,6 @@ class Chip(FocalPlane): self.fdModel = None else: try: - # with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion: with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion: with open(field_distortion, "rb") as f: self.fdModel = pickle.load(f) @@ -88,9 +83,11 @@ class Chip(FocalPlane): # Get boundary (in pix) self.bound = self.getChipLim() + self.ccdEffCurve_dir = ccdEffCurve_dir self.CRdata_dir = CRdata_dir - slsconfs = self.getChipSLSConf() + + slsconfs = chip_utils.getChipSLSConf(chipID=self.chipID) if np.size(slsconfs) == 1: try: with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs) as conf_path: @@ -157,13 +154,6 @@ class Chip(FocalPlane): if filter_type in ['NUV', 'u', 'GU']: filename = 'UV0.txt' if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt' - # Mirror efficiency: - # if filter_type == 'NUV': mirror_eff = 0.54 - # if filter_type == 'u': mirror_eff = 0.68 - # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8 - # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right - # path = os.path.join(self.ccdEffCurve_dir, filename) - # table = Table.read(path, format='ascii') try: with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filename) as ccd_path: table = Table.read(ccd_path, format='ascii') @@ -181,13 +171,26 @@ class Chip(FocalPlane): except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path: self.attachedSizes = np.loadtxt(cr_path) + + def loadSLSFLATCUBE(self, flat_fn='flat_cube.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 - def getChipFilter(self, chipID=None, filter_layout=None): + def getChipFilter(self, chipID=None): """Return the filter index and type for a given chip #(chipID) """ filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"] - if filter_layout is not None: - return filter_layout[chipID][0], filter_layout[chipID][1] if chipID == None: chipID = self.chipID @@ -217,45 +220,18 @@ class Chip(FocalPlane): Returns: A galsim BoundsD object """ - if ((chipID is not None) and (int(chipID) <= 30)) or (self.chipID <= 30): - # [TODO] - if chipID == None: - chipID = self.chipID - rowID, colID = self.rowID, self.colID - else: - rowID, colID = self.getChipRowCol(chipID) - gx1, gx2 = self.npix_gap_x - gy = self.npix_gap_y - - # xlim of a given CCD chip - xrem = 2*(colID - 1) - (self.nchip_x - 1) - xcen = (self.npix_x//2 + gx1//2) * xrem - if chipID >= 26 or chipID == 21: - xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1) - if chipID <= 5 or chipID == 10: - xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1) - nx0 = xcen - self.npix_x//2 + 1 - nx1 = xcen + self.npix_x//2 - - # ylim of a given CCD chip - yrem = (rowID - 1) - self.nchip_y // 2 - ycen = (self.npix_y + gy) * yrem - ny0 = ycen - self.npix_y//2 + 1 - ny1 = ycen + self.npix_y//2 - return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1) - else: - xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10 - xcen = self.x_cen / self.pix_size - ycen = self.y_cen / self.pix_size - sign_x = [-1., 1., -1., 1.] - sign_y = [-1., -1., 1., 1.] - for i in range(4): - x = xcen + sign_x[i] * self.npix_x / 2. - y = ycen + sign_y[i] * self.npix_y / 2. - x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle) - xmin, xmax = min(xmin, x), max(xmax, x) - ymin, ymax = min(ymin, y), max(ymax, y) - return galsim.BoundsD(xmin, xmax, ymin, ymax) + 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) @@ -291,83 +267,8 @@ class Chip(FocalPlane): noise = self.dark_noise * exptime + self.read_noise**2 return noise - def getChipSLSConf(self): - confFile = '' - if self.chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf'] - if self.chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf'] - if self.chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf'] - if self.chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf'] - if self.chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf'] - if self.chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf'] - if self.chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf'] - if self.chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf'] - if self.chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf'] - if self.chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf'] - if self.chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf'] - if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf'] - return confFile - - def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200): - datetime_obs = datetime.utcfromtimestamp(timestamp) - date_obs = datetime_obs.strftime("%y%m%d") - time_obs = datetime_obs.strftime("%H%M%S") - h_prim = generatePrimaryHeader( - xlen=self.npix_x, - ylen=self.npix_y, - pointNum = str(pointing_ID), - ra=ra_cen, - dec=dec_cen, - pixel_scale=self.pix_scale, - date=date_obs, - time_obs=time_obs, - im_type = im_type, - exptime=exptime, - chip_name=str(self.chipID).rjust(2, '0') - ) - h_ext = generateExtensionHeader( - chip=self, - xlen=self.npix_x, - ylen=self.npix_y, - ra=ra_cen, - dec=dec_cen, - pa=img_rot.deg, - gain=self.gain, - readout=self.read_noise, - dark=self.dark_noise, - saturation=90000, - pixel_scale=self.pix_scale, - pixel_size=self.pix_size, - xcen=self.x_cen, - ycen=self.y_cen, - extName='SCI', - timestamp = timestamp, - exptime = exptime, - readoutTime = 40.) - return h_prim, h_ext - - def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200): - 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, - exptime=exptime, - timestamp = timestamp) - hdu1 = fits.PrimaryHDU(header=h_prim) - hdu1.add_checksum() - hdu1.header.comments['CHECKSUM'] = 'HDU checksum' - hdu1.header.comments['DATASUM'] = 'data unit checksum' - hdu2 = fits.ImageHDU(img.array, header=h_ext) - hdu2.add_checksum() - hdu2.header.comments['XTENSION'] = 'extension type' - hdu2.header.comments['CHECKSUM'] = 'HDU checksum' - hdu2.header.comments['DATASUM'] = 'data unit checksum' - hdu1 = fits.HDUList([hdu1, hdu2]) - fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits') - hdu1.writeto(fname, output_verify='ignore', overwrite=True) - def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None): + # Set random seeds SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"]) @@ -386,42 +287,19 @@ class Chip(FocalPlane): self.logger = logger # Get Poisson noise generator - seed = int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID - rng_poisson = galsim.BaseDeviate(seed) - poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.) + rng_poisson, poisson_noise = chip_utils.get_poisson( + seed=int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID, sky_level=0.) # Add sky background - if sky_map is None: - sky_map = filt.getSkyNoise(exptime=exptime) - sky_map = sky_map * np.ones_like(img.array) - sky_map = galsim.Image(array=sky_map) - # Apply Poisson noise to the sky map - # (NOTE): only for photometric chips - # since it utilize the photon shooting - # to draw stamps - if self.survey_type == "photometric": - sky_map.addNoise(poisson_noise) - elif img.array.shape != sky_map.shape: - raise ValueError("The shape img and sky_map must be equal.") - elif tel is not None: # If sky_map is given in flux - sky_map = sky_map * tel.pupil_area * exptime if config["ins_effects"]["add_back"] == True: - img += sky_map - del sky_map + img, sky_map = chip_utils.add_sky_background(img=img, filt=filt, exptime=exptime, sky_map=sky_map, tel=tel) + del sky_map # Apply flat-field large scale structure for one chip if config["ins_effects"]["flat_fielding"] == True: - if self.logger is not None: - self.logger.info(" Creating and applying Flat-Fielding") - msg = str(img.bounds) - self.logger.info(msg) - else: - print(" Creating and applying Flat-Fielding", flush=True) - print(img.bounds, flush=True) - flat_img = effects.MakeFlatSmooth( - img.bounds, - int(config["random_seeds"]["seed_flat"])) - flat_normal = flat_img / np.mean(flat_img.array) + chip_utils.log_info(msg=" Creating and applying Flat-Fielding", logger=self.logger) + chip_utils.log_info(msg=str(img.bounds), logger=self.logger) + flat_img, flat_normal = chip_utils.get_flat(img=img, seed=int(config["random_seeds"]["seed_flat"])) if self.survey_type == "photometric": img *= flat_normal del flat_normal @@ -430,10 +308,7 @@ class Chip(FocalPlane): # Apply Shutter-effect for one chip if config["ins_effects"]["shutter_effect"] == True: - if self.logger is not None: - self.logger.info(" Apply shutter effect") - else: - print(" Apply shutter effect", flush=True) + chip_utils.log_info(msg=" Apply shutter effect", logger=self.logger) shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3) # shutter effect normalized image for this chip if self.survey_type == "photometric": img *= shuttimg @@ -443,36 +318,19 @@ class Chip(FocalPlane): del shutt_gsimg del shuttimg - # Add Poisson noise to the resulting images - # (NOTE): this can only applied to the slitless image - # since it dose not use photon shooting to draw stamps - if self.survey_type == "spectroscopic": - img.addNoise(poisson_noise) + # # Add Poisson noise to the resulting images + # # (NOTE): this can only applied to the slitless image + # # since it dose not use photon shooting to draw stamps + # if self.survey_type == "spectroscopic": + # img.addNoise(poisson_noise) # Add cosmic-rays if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS': - if self.logger is not None: - self.logger.info((" Adding Cosmic-Ray")) - else: - print(" Adding Cosmic-Ray", flush=True) - cr_map, cr_event_num = effects.produceCR_Map( - xLen=self.npix_x, yLen=self.npix_y, - exTime=exptime+0.5*self.readout_time, - cr_pixelRatio=0.003*(exptime+0.5*self.readout_time)/600., - gain=self.gain, - attachedSizes=self.attachedSizes, - seed=SeedCosmicRay+pointing_ID*30+self.chipID) # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3; - img += cr_map - cr_map[cr_map > 65535] = 65535 - cr_map[cr_map < 0] = 0 - crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16) - del cr_map - # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID)) - # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID)) - # datetime_obs = datetime.utcfromtimestamp(timestamp_obs) - # date_obs = datetime_obs.strftime("%y%m%d") - # time_obs = datetime_obs.strftime("%H%M%S") - self.outputCal( + chip_utils.log_info(msg=" Adding Cosmic-Ray", logger=self.logger) + img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, + seed=SeedCosmicRay+pointing_ID*30+self.chipID) + chip_utils.outputCal( + chip=self, img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, @@ -486,25 +344,28 @@ class Chip(FocalPlane): # Apply PRNU effect and output PRNU flat file: if config["ins_effects"]["prnu_effect"] == True: - if self.logger is not None: - self.logger.info(" Applying PRNU effect") - else: - print(" Applying PRNU effect", flush=True) - prnu_img = effects.PRNU_Img( - xsize=self.npix_x, - ysize=self.npix_y, - sigma=0.01, - seed=int(config["random_seeds"]["seed_prnu"]+self.chipID)) - img *= prnu_img + chip_utils.log_info(msg=" Applying PRNU effect", logger=self.logger) + img, prnu_img = chip_utils.add_PRNU(img=img, chip=self, + seed=int(config["random_seeds"]["seed_prnu"]+self.chipID)) if config["output_setting"]["prnu_output"] == True: prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID)) if config["output_setting"]["flat_output"] == False: del prnu_img - # Add dark current + # # Add dark current + # if config["ins_effects"]["add_dark"] == True: + # dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time))) + # img.addNoise(dark_noise) + + # Add dark current & Poisson noise if config["ins_effects"]["add_dark"] == True: - dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time))) - img.addNoise(dark_noise) + img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise) + else: + img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise, dark_noise=0.) + + # Add diffusion & brighter-fatter effects + if config["ins_effects"]["bright_fatter"] == True: + img = chip_utils.add_brighter_fatter(img=img) # Add Hot Pixels or/and Dead Pixels rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID))) @@ -517,34 +378,22 @@ class Chip(FocalPlane): # Apply Nonlinearity on the chip image if config["ins_effects"]["non_linear"] == True: - if self.logger is not None: - self.logger.info(" Applying Non-Linearity on the chip image") - else: - print(" Applying Non-Linearity on the chip image", flush=True) + chip_utils.log_info(msg=" Applying Non-Linearity on the chip image", logger=self.logger) img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0) # Apply CCD Saturation & Blooming if config["ins_effects"]["saturbloom"] == True: - if self.logger is not None: - self.logger.info(" Applying CCD Saturation & Blooming") - else: - print(" Applying CCD Saturation & Blooming") + chip_utils.log_info(msg=" Applying CCD Saturation & Blooming", logger=self.logger) img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell) # Apply CTE Effect if config["ins_effects"]["cte_trail"] == True: - if self.logger is not None: - self.logger.info(" Apply CTE Effect") - else: - print(" Apply CTE Effect") + chip_utils.log_info(msg=" Apply CTE Effect", logger=self.logger) img = effects.CTE_Effect(GSImage=img, threshold=27) # Add Bias level if config["ins_effects"]["add_bias"] == True: - if self.logger is not None: - self.logger.info(" Adding Bias level and 16-channel non-uniformity") - else: - print(" Adding Bias level and 16-channel non-uniformity") + chip_utils.log_info(msg=" Adding Bias level and 16-channel non-uniformity", logger=self.logger) if config["ins_effects"]["bias_16channel"] == True: img = effects.AddBiasNonUniform16(img, bias_level=float(self.bias_level), @@ -562,10 +411,7 @@ class Chip(FocalPlane): img.addNoise(readout_noise) # Apply Gain & Quantization - if self.logger is not None: - self.logger.info(" Applying Gain (and 16 channel non-uniformity) & Quantization") - else: - print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True) + chip_utils.log_info(msg=" Applying Gain (and 16 channel non-uniformity) & Quantization", logger=self.logger) if config["ins_effects"]["gain_16channel"] == True: img, self.gain_channel = effects.ApplyGainNonUniform16( img, gain=self.gain, @@ -633,12 +479,9 @@ class Chip(FocalPlane): BiasCombImg.replaceNegative(replace_value=0) BiasCombImg.quantize() BiasCombImg = galsim.ImageUS(BiasCombImg) - # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1)) - # datetime_obs = datetime.utcfromtimestamp(timestamp_obs) - # date_obs = datetime_obs.strftime("%y%m%d") - # time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 10 * 60 - self.outputCal( + chip_utils.outputCal( + chip=self, img=BiasCombImg, ra_cen=ra_cen, dec_cen=dec_cen, @@ -736,12 +579,9 @@ class Chip(FocalPlane): FlatCombImg.replaceNegative(replace_value=0) FlatCombImg.quantize() FlatCombImg = galsim.ImageUS(FlatCombImg) - # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1)) - # datetime_obs = datetime.utcfromtimestamp(timestamp_obs) - # date_obs = datetime_obs.strftime("%y%m%d") - # time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 10 * 60 - self.outputCal( + chip_utils.outputCal( + chip=self, img=FlatCombImg, ra_cen=ra_cen, dec_cen=dec_cen, @@ -793,10 +633,8 @@ class Chip(FocalPlane): cr_map[cr_map < 0] = 0 crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16) del cr_map - # datetime_obs = datetime.utcfromtimestamp(timestamp_obs) - # date_obs = datetime_obs.strftime("%y%m%d") - # time_obs = datetime_obs.strftime("%H%M%S") - self.outputCal( + chip_utils.outputCal( + chip=self, img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, @@ -860,12 +698,9 @@ class Chip(FocalPlane): DarkCombImg.replaceNegative(replace_value=0) DarkCombImg.quantize() DarkCombImg = galsim.ImageUS(DarkCombImg) - # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1)) - # datetime_obs = datetime.utcfromtimestamp(timestamp_obs) - # date_obs = datetime_obs.strftime("%y%m%d") - # time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 10 * 60 - self.outputCal( + chip_utils.outputCal( + chip=chip, img=DarkCombImg, ra_cen=ra_cen, dec_cen=dec_cen, @@ -894,19 +729,3 @@ class Chip(FocalPlane): # del sub_img return img - def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'): - from astropy.io import fits - try: - with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path: - flat_fits = fits.open(data_path, ignore_missing_simple=True) - except AttributeError: - with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path: - flat_fits = fits.open(data_path, ignore_missing_simple=True) - - fl = len(flat_fits) - fl_sh = flat_fits[0].data.shape - assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!' - self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]]) - for i in np.arange(0, fl, 1): - self.flat_cube[i] = flat_fits[i].data - diff --git a/ObservationSim/Instrument/Chip/ChipUtils.py b/ObservationSim/Instrument/Chip/ChipUtils.py new file mode 100644 index 0000000..97825a7 --- /dev/null +++ b/ObservationSim/Instrument/Chip/ChipUtils.py @@ -0,0 +1,196 @@ +import os +import galsim +import ctypes +import numpy as np +from astropy.io import fits +from datetime import datetime + +from ObservationSim.Instrument.Chip import Effects as effects +from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader + +try: + import importlib.resources as pkg_resources +except ImportError: + # Try backported to PY<37 'importlib_resources' + import importlib_resources as pkg_resources + + +def log_info(msg, logger=None): + if logger: + logger.info(msg) + else: + print(msg, flush=True) + +def getChipSLSConf(chipID): + confFile = '' + if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf'] + if chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf'] + if chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf'] + if chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf'] + if chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf'] + if chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf'] + if chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf'] + if chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf'] + if chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf'] + if chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf'] + if chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf'] + if chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf'] + return confFile + +def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200): + datetime_obs = datetime.utcfromtimestamp(timestamp) + date_obs = datetime_obs.strftime("%y%m%d") + time_obs = datetime_obs.strftime("%H%M%S") + h_prim = generatePrimaryHeader( + xlen=chip.npix_x, + ylen=chip.npix_y, + pointNum = str(pointing_ID), + ra=ra_cen, + dec=dec_cen, + pixel_scale=chip.pix_scale, + date=date_obs, + time_obs=time_obs, + im_type = im_type, + exptime=exptime, + chip_name=str(chip.chipID).rjust(2, '0') + ) + h_ext = generateExtensionHeader( + chip=chip, + xlen=chip.npix_x, + ylen=chip.npix_y, + ra=ra_cen, + dec=dec_cen, + pa=img_rot.deg, + gain=chip.gain, + readout=chip.read_noise, + dark=chip.dark_noise, + saturation=90000, + pixel_scale=chip.pix_scale, + pixel_size=chip.pix_size, + xcen=chip.x_cen, + ycen=chip.y_cen, + extName='SCI', + timestamp = timestamp, + exptime = exptime, + readoutTime = chip.readout_time) + return h_prim, h_ext + +def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200): + h_prim, h_ext = generateHeader( + chip=chip, + ra_cen=ra_cen, + dec_cen=dec_cen, + img_rot=img_rot, + im_type=im_type, + pointing_ID=pointing_ID, + exptime=exptime, + timestamp=timestamp) + hdu1 = fits.PrimaryHDU(header=h_prim) + hdu1.add_checksum() + hdu1.header.comments['CHECKSUM'] = 'HDU checksum' + hdu1.header.comments['DATASUM'] = 'data unit checksum' + hdu2 = fits.ImageHDU(img.array, header=h_ext) + hdu2.add_checksum() + hdu2.header.comments['XTENSION'] = 'extension type' + hdu2.header.comments['CHECKSUM'] = 'HDU checksum' + hdu2.header.comments['DATASUM'] = 'data unit checksum' + hdu1 = fits.HDUList([hdu1, hdu2]) + fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits') + hdu1.writeto(fname, output_verify='ignore', overwrite=True) + +def add_sky_background(img, filt, exptime, sky_map=None, tel=None): + # Add sky background + if sky_map is None: + sky_map = filt.getSkyNoise(exptime=exptime) + sky_map = sky_map * np.ones_like(img.array) + sky_map = galsim.Image(array=sky_map) + # Apply Poisson noise to the sky map + # # (NOTE): only for photometric chips if it utilizes the photon shooting to draw stamps + # if self.survey_type == "photometric": + # sky_map.addNoise(poisson_noise) + elif img.array.shape != sky_map.shape: + raise ValueError("The shape img and sky_map must be equal.") + elif tel is not None: # If sky_map is given in flux + sky_map = sky_map * tel.pupil_area * exptime + img += sky_map + return img, sky_map + +def get_flat(img, seed): + flat_img = effects.MakeFlatSmooth( + GSBounds=img.bounds, + seed=seed) + flat_normal = flat_img / np.mean(flat_img.array) + return flat_img, flat_normal + +def add_cosmic_rays(img, chip, exptime=150, seed=0): + cr_map, cr_event_num = effects.produceCR_Map( + xLen=chip.npix_x, yLen=chip.npix_y, + exTime=exptime+0.5*chip.readout_time, + cr_pixelRatio=0.003*(exptime+0.5*chip.readout_time)/600., + gain=chip.gain, + attachedSizes=chip.attachedSizes, + seed=seed) # 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 + return img, crmap_gsimg, cr_event_num + +def add_PRNU(img, chip, seed=0): + prnu_img = effects.PRNU_Img( + xsize=chip.npix_x, + ysize=chip.npix_y, + sigma=0.01, + seed=seed) + img *= prnu_img + return img, prnu_img + +def get_poisson(seed=0, sky_level=0.): + rng_poisson = galsim.BaseDeviate(seed) + poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=sky_level) + return rng_poisson, poisson_noise + +def get_base_img(img, read_noise, readout_time, dark_noise, exptime=150.): + base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time) + base_img = base_level * np.ones_like(img.array) + return base_img + +def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=None, dark_noise=None): + if poisson_noise is None: + _, poisson_noise = get_poisson(seed=seed, sky_level=sky_level) + read_noise = chip.read_noise + if dark_noise is None: + dark_noise = chip.dark_noise + base_img = get_base_img(img=img, read_noise=read_noise, readout_time=chip.readout_time, dark_noise=dark_noise, exptime=exptime) + img += base_img + img.addNoise(poisson_noise) + img -= read_noise**2 + return img, base_img + +def add_brighter_fatter(img): + #Inital dynamic lib + try: + with pkg_resources.files('ObservationSim.Instrument.Chip.lib_bf').joinpath("libmoduleBF.so") as lib_path: + lib_bf = ctypes.CDLL(lib_path) + except AttributeError: + with pkg_resources.path('ObservationSim.Instrument.Chip.lib_bf', "libmoduleBF.so") as lib_path: + lib_bf = ctypes.CDLL(lib_path) + lib_bf.addEffects.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_int] + + # Set bit flag + bit_flag = 1 + bit_flag = bit_flag | (1 << 2) + + nx, ny = img.array.shape + nn = nx * ny + arr_ima= (ctypes.c_float*nn)() + arr_imc= (ctypes.c_float*nn)() + + arr_ima[:]= img.array.reshape(nn) + arr_imc[:]= np.zeros(nn) + + lib_bf.addEffects(nx, ny, arr_ima, arr_imc, bit_flag) + img.array[:, :] = np.array(arr_imc[:]).reshape([nx, ny]) + del arr_ima, arr_imc + return img \ No newline at end of file diff --git a/ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c b/ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c new file mode 100644 index 0000000..199b48c --- /dev/null +++ b/ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c @@ -0,0 +1,120 @@ +#include +#include +#include +#include +#include "nrutil.h" + +#define ISSETBITFLAG(x,b) ((x) & (1 << (b))) +#define ADD_DIFFUSION 1 +#define ADD_BF_FILTER 2 + +float linearInterp(float xp, float x0, float y0, float x1, float y1); + +void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bit_flag) +{ + int nx, ny, i,j,k,ks; + int it,jt,itt,jtt; + int diffuidx[26][2],diffuN,ilow,ih,im,dim[3]; + float diffua[5][5],cdiffu[26],**bfa; + double mvar,mcov,tmp,ma,mb,mc; + char fname[100]; + + nx = ngx_ima; //input-image size + ny = ngy_ima; + + //0. init. original image with an input array (arr_ima) + //1. Adding diffusion effect. + if(ISSETBITFLAG(bit_flag, ADD_DIFFUSION)) + { + printf("adding diffusion.....\n"); + printf("ERR: no diffusion filter ..."); + exit(0); + } + + + //2. Adding BF effect + if(ISSETBITFLAG(bit_flag, ADD_BF_FILTER)) + { + printf("Adding BF effect...\n"); + //setup BF correlation fliter + float neX; + float neP1 = 50000; + float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302}; + float neP2 = 10000; + float bfaP2[9]={0.9945288, 0.0003041936, 0.0007539311, 0.0002424675, 0.0001226098, 0.00009308617, 0.00008027447, 0.00006309676, 0.00006400052}; + + bfa=matrix(-2,2,-2,2); + + // smooth with the BF filter + for(i=0;i= 10000) + { + bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0; + bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575; + bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652; + bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335; + bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]); + bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118; + bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]); + bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083; + bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043; + } + else + { + neX=10000; + bfa[0][0]=0; + bfa[0][1]=bfa[0][-1]=bfaP2[1]; + bfa[-1][0]=bfa[1][0]=bfaP2[2]; + bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=bfaP2[3]; + bfa[0][-2]=bfa[0][2]=bfaP2[4]; + bfa[-2][0]=bfa[2][0]=bfaP2[5]; + bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=bfaP2[6]; + bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=bfaP2[7]; + bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=bfaP2[8]; + } + tmp = 0; + for(it=-2;it<=2;it++) + for(jt=-2;jt<=2;jt++) + { + bfa[it][jt] = bfa[it][jt]/neX*arr_ima[j+i*ny]; + tmp += bfa[it][jt]; + } + bfa[0][0]=1.-tmp; + + // assign electrons according to the BF filter bfat + for(it=-2;it<=2;it++) + { + for(jt=-2;jt<=2;jt++) + { + itt=i+it; + jtt=j+jt; + if(itt>=0 && jtt>=0 && itt +#include +#include +#define NR_END 1 +#define FREE_ARG char* + +void nrerror(char error_text[]) +/* Numerical Recipes standard error handler */ +{ + fprintf(stderr,"Numerical Recipes run-time error...\n"); + fprintf(stderr,"%s\n",error_text); + fprintf(stderr,"...now exiting to system...\n"); + exit(1); +} + +float *vector(long nl, long nh) +/* allocate a float vector with subscript range v[nl..nh] */ +{ + float *v; + + v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); + if (!v) nrerror("allocation failure in vector()"); + return v-nl+NR_END; +} + +int *ivector(long nl, long nh) +/* allocate an int vector with subscript range v[nl..nh] */ +{ + int *v; + + v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); + if (!v) nrerror("allocation failure in ivector()"); + return v-nl+NR_END; +} + +unsigned char *cvector(long nl, long nh) +/* allocate an unsigned char vector with subscript range v[nl..nh] */ +{ + unsigned char *v; + + v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); + if (!v) nrerror("allocation failure in cvector()"); + return v-nl+NR_END; +} + +long *lvector(long nl, long nh) +/* allocate an long vector with subscript range v[nl..nh] */ +{ + long *v; + + v=(long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); + if (!v) nrerror("allocation failure in lvector()"); + return v-nl+NR_END; +} + +double *dvector(long nl, long nh) +/* allocate a double vector with subscript range v[nl..nh] */ +{ + double *v; + + v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); + if (!v) nrerror("allocation failure in dvector()"); + return v-nl+NR_END; +} + +float **matrix(long nrl, long nrh, long ncl, long nch) +/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + float **m; + + /* allocate pointers to rows */ + m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +double **dmatrix(long nrl, long nrh, long ncl, long nch) +/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + double **m; + + /* allocate pointers to rows */ + m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +int **imatrix(long nrl, long nrh, long ncl, long nch) +/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + int **m; + + /* allocate pointers to rows */ + m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + + /* allocate rows and set pointers to them */ + m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, + long newrl, long newcl) +/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ +{ + long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; + float **m; + + /* allocate array of pointers to rows */ + m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure in submatrix()"); + m += NR_END; + m -= newrl; + + /* set pointers to rows */ + for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch) +/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix +declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 +and ncol=nch-ncl+1. The routine should be called with the address +&a[0][0] as the first argument. */ +{ + long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; + float **m; + + /* allocate pointers to rows */ + m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure in convert_matrix()"); + m += NR_END; + m -= nrl; + + /* set pointers to rows */ + m[nrl]=a-ncl; + for(i=1,j=nrl+1;i +#define NR_END 1 +#define FREE_ARG char* + +void nrerror(error_text) +char error_text[]; +/* Numerical Recipes standard error handler */ +{ + void exit(); + + fprintf(stderr,"Numerical Recipes run-time error...\n"); + fprintf(stderr,"%s\n",error_text); + fprintf(stderr,"...now exiting to system...\n"); + exit(1); +} + +float *vector(nl,nh) +long nh,nl; +/* allocate a float vector with subscript range v[nl..nh] */ +{ + float *v; + + v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float))); + if (!v) nrerror("allocation failure in vector()"); + return v-nl+NR_END; +} + +int *ivector(nl,nh) +long nh,nl; +/* allocate an int vector with subscript range v[nl..nh] */ +{ + int *v; + + v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int))); + if (!v) nrerror("allocation failure in ivector()"); + return v-nl+NR_END; +} + +unsigned char *cvector(nl,nh) +long nh,nl; +/* allocate an unsigned char vector with subscript range v[nl..nh] */ +{ + unsigned char *v; + + v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char))); + if (!v) nrerror("allocation failure in cvector()"); + return v-nl+NR_END; +} + +long *lvector(nl,nh) +long nh,nl; +/* allocate an unsigned long vector with subscript range v[nl..nh] */ +{ + long *v; + + v=(long *)malloc((int) ((nh-nl+1+NR_END)*sizeof(long))); + if (!v) nrerror("allocation failure in lvector()"); + return v-nl+NR_END; +} + +double *dvector(nl,nh) +long nh,nl; +/* allocate a double vector with subscript range v[nl..nh] */ +{ + double *v; + + v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double))); + if (!v) nrerror("allocation failure in dvector()"); + return v-nl+NR_END; +} + +float **matrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + float **m; + + /* allocate pointers to rows */ + m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +double **dmatrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + double **m; + + /* allocate pointers to rows */ + m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +int **imatrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + int **m; + + /* allocate pointers to rows */ + m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + + /* allocate rows and set pointers to them */ + m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl) +float **a; +long newcl,newrl,oldch,oldcl,oldrh,oldrl; +/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ +{ + long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; + float **m; + + /* allocate array of pointers to rows */ + m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure in submatrix()"); + m += NR_END; + m -= newrl; + + /* set pointers to rows */ + for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +float **convert_matrix(a,nrl,nrh,ncl,nch) +float *a; +long nch,ncl,nrh,nrl; +/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix +declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 +and ncol=nch-ncl+1. The routine should be called with the address +&a[0][0] as the first argument. */ +{ + long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; + float **m; + + /* allocate pointers to rows */ + m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("allocation failure in convert_matrix()"); + m += NR_END; + m -= nrl; + + /* set pointers to rows */ + m[nrl]=a-ncl; + for(i=1,j=nrl+1;i (dmaxarg2) ?\ + (dmaxarg1) : (dmaxarg2)) + +static double dminarg1,dminarg2; +#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\ + (dminarg1) : (dminarg2)) + +static float maxarg1,maxarg2; +#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\ + (maxarg1) : (maxarg2)) + +static float minarg1,minarg2; +#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\ + (minarg1) : (minarg2)) + +static long lmaxarg1,lmaxarg2; +#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\ + (lmaxarg1) : (lmaxarg2)) + +static long lminarg1,lminarg2; +#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\ + (lminarg1) : (lminarg2)) + +static int imaxarg1,imaxarg2; +#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\ + (imaxarg1) : (imaxarg2)) + +static int iminarg1,iminarg2; +#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\ + (iminarg1) : (iminarg2)) + +#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) + +void nrerror(char error_text[]); +float *vector(long nl, long nh); +int *ivector(long nl, long nh); +unsigned char *cvector(long nl, long nh); +long *lvector(long nl, long nh); +double *dvector(long nl, long nh); +float **matrix(long nrl, long nrh, long ncl, long nch); +double **dmatrix(long nrl, long nrh, long ncl, long nch); +int **imatrix(long nrl, long nrh, long ncl, long nch); +float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, + long newrl, long newcl); +float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch); +float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh); +void free_vector(float *v, long nl, long nh); +void free_ivector(int *v, long nl, long nh); +void free_cvector(unsigned char *v, long nl, long nh); +void free_lvector(long *v, long nl, long nh); +void free_dvector(double *v, long nl, long nh); +void free_matrix(float **m, long nrl, long nrh, long ncl, long nch); +void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch); +void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch); +void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch); +void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch); +void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch, + long ndl, long ndh); + + +int ***i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh); +void free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch, + long ndl, long ndh); + +unsigned char ***b3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh); +void free_b3tensor(unsigned char ***t, long nrl, long nrh, long ncl, long nch, + long ndl, long ndh); + +double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh); +void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch, + long ndl, long ndh); + + +char **cmatrix(long nrl, long nrh, long ncl, long nch); +void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch); + +#endif /* _NR_UTILS_H_ */ diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 6bda701..81ce517 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -14,18 +14,6 @@ from ObservationSim.MockObject.MockObject import MockObject class Galaxy(MockObject): def __init__(self, param, logger=None): super().__init__(param, logger=logger) - # self.thetaR = self.param["theta"] - # self.bfrac = self.param["bfrac"] - # self.hlr_disk = self.param["hlr_disk"] - # self.hlr_bulge = self.param["hlr_bulge"] - - # Extract ellipticity components - # self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees) - # self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees) - # self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees) - # self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2 - # self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2 - # self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2 if not hasattr(self, "disk_sersic_idx"): self.disk_sersic_idx = 1. @@ -101,10 +89,6 @@ class Galaxy(MockObject): if self.logger: self.logger.error(e) return 2, None - - nphotons_sum = 0 - photons_list = [] - xmax, ymax = 0, 0 # # [C6 TEST] # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge)) @@ -114,7 +98,7 @@ class Galaxy(MockObject): if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy big_galaxy = True - # (TEST) Galsim Parameters + # Set Galsim Parameters if self.getMagFilter(filt) <= 15 and (not big_galaxy): folding_threshold = 5.e-4 else: @@ -122,7 +106,7 @@ class Galaxy(MockObject): gsp = galsim.GSParams(folding_threshold=folding_threshold) self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, - img_real_wcs=self.real_wcs) + img_real_wcs=self.chip_wcs) x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) @@ -130,9 +114,11 @@ class Galaxy(MockObject): dx = x - x_nominal dy = y - y_nominal offset = galsim.PositionD(dx, dy) + # Get real local wcs of object (deal with chip rotation w.r.t its center) + chip_wcs_local = self.chip_wcs.local(self.real_pos) + is_updated = 0 - real_wcs_local = self.real_wcs.local(self.real_pos) - + # Model the galaxy as disk + bulge disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp) disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) disk = disk.shear(disk_shape) @@ -140,33 +126,33 @@ class Galaxy(MockObject): bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge = bulge.shear(bulge_shape) + # Get shear and deal with shear induced by field distortion if fd_shear: g1 += fd_shear.g1 g2 += fd_shear.g2 gal_shear = galsim.Shear(g1=g1, g2=g2) + # Loop over all sub-bandpasses for i in range(len(bandpass_list)): bandpass = bandpass_list[i] - try: sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) if self.logger: self.logger.error(e) - # return False continue - ratio = sub/full if not (ratio == -1 or (ratio != ratio)): nphotons = ratio * nphotons_tot else: continue - nphotons_sum += nphotons - + + # nphotons_sum += nphotons # # [C6 TEST] # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) + # Get PSF model psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk @@ -194,39 +180,21 @@ class Galaxy(MockObject): # for stat in top_stats[:10]: # print(stat) - stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) - photons = stamp.photons - photons.x += x_nominal - photons.y += y_nominal - photons_list.append(photons) - - stamp.wcs = real_wcs_local + # stamp = gal.drawImage(wcs=chip_wcs_local, method='phot', offset=offset, save_photons=True) + stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset) + if np.sum(np.isnan(stamp.array)) > 0: + # ERROR happens + return 2, pos_shear stamp.setCenter(x_nominal, y_nominal) bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) - if bounds.area() > 0: chip.img.setOrigin(0, 0) - stamp[bounds] = chip.img[bounds] - - if not big_galaxy: - for i in range(len(photons_list)): - if i == 0: - chip.sensor.accumulate(photons_list[i], stamp) - else: - chip.sensor.accumulate(photons_list[i], stamp, resume=True) - else: - sensor = galsim.Sensor() - for i in range(len(photons_list)): - if i == 0: - sensor.accumulate(photons_list[i], stamp) - else: - sensor.accumulate(photons_list[i], stamp, resume=True) - del sensor - - chip.img[bounds] = stamp[bounds] - + chip.img[bounds] += stamp[bounds] + is_updated = 1 chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) - else: + del stamp + + if is_updated == 0: # Return code 0: object photons missed this detector print("obj %s missed"%(self.id)) if self.logger: @@ -235,8 +203,6 @@ class Galaxy(MockObject): # # [C6 TEST] # print("nphotons_sum = ", nphotons_sum) - del photons_list - del stamp return 1, pos_shear def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, @@ -255,7 +221,7 @@ class Galaxy(MockObject): names=('WAVELENGTH', 'FLUX')) self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, - img_real_wcs=self.real_wcs) + img_real_wcs=self.chip_wcs) x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) @@ -264,7 +230,7 @@ class Galaxy(MockObject): dy = y - y_nominal offset = galsim.PositionD(dx, dy) - real_wcs_local = self.real_wcs.local(self.real_pos) + chip_wcs_local = self.chip_wcs.local(self.real_pos) big_galaxy = False @@ -313,7 +279,7 @@ class Galaxy(MockObject): # if fd_shear is not None: # gal = gal.shear(fd_shear) - starImg = gal.drawImage(wcs=real_wcs_local, offset=offset) + starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset) origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_nominal - (starImg.center.x - starImg.xmin)] @@ -340,7 +306,7 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) @@ -357,7 +323,7 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) del sdp_p1 del sdp_p2 @@ -369,7 +335,7 @@ class Galaxy(MockObject): conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) del sdp elif grating_split_pos_chip>=gal_end[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, @@ -379,7 +345,7 @@ class Galaxy(MockObject): conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) @@ -404,46 +370,6 @@ class Galaxy(MockObject): final = galsim.Convolve(psf, gal) return final - def drawObject(self, img, final, noise_level=0.0, flux=None, filt=None, tel=None, exptime=150.): - """ Override the method in parent class - Need to constrain the size of image stamp for extended objects - """ - isUpdated = True - if flux == None: - flux = self.getElectronFluxFilt(filt, tel, exptime) - stamp = final.drawImage(wcs=self.localWCS, offset=self.offset) - stamp_arr = stamp.array - mask = (stamp_arr >= 0.001*noise_level) # why 0.001? - err = int(np.sqrt(mask.sum())) - if np.mod(err, 2) == 1: - err += 1 - # if err == 1: - if err == 0: - subSize = 16 # why 16? - else: - subSize = max([err, 16]) - fluxRatio = flux / stamp_arr[mask].sum() - final = final.withScaledFlux(fluxRatio) - - imgSub = galsim.ImageF(subSize, subSize) - - # Draw with FFT - # stamp = final.drawImage(image=imgSub, wcs=self.localWCS, offset=self.offset) - - # Draw with Photon Shoot - stamp = final.drawImage(image=imgSub, wcs=self.localWCS, method='phot', offset=self.offset) - - stamp.setCenter(self.x_nominal, self.y_nominal) - if np.sum(np.isnan(stamp.array)) >= 1: - stamp.setZero() - bounds = stamp.bounds & img.bounds - if bounds.area() == 0: - isUpdated = False - else: - img[bounds] += stamp[bounds] - return img, stamp, isUpdated - - def getObservedEll(self, g1=0, g2=0): e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2) return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index ae7b708..654cb72 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -13,7 +13,6 @@ from ObservationSim.MockObject.SpecDisperser import SpecDisperser class MockObject(object): def __init__(self, param, logger=None): self.param = param - for key in self.param: setattr(self, key, self.param[key]) @@ -27,14 +26,10 @@ class MockObject(object): elif self.param["star"] == 3: self.type = "stamp" ###mock_stamp_END - self.sed = None - + self.fd_shear = None # Place holder for outputs self.additional_output_str = "" - - self.fd_shear = None - self.logger = logger def getMagFilter(self, filt): @@ -65,6 +60,7 @@ class MockObject(object): def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None): self.posImg = img.wcs.toImage(self.getPosWorld()) self.localWCS = img.wcs.local(self.posImg) + # Apply field distortion model if (fdmodel is not None) and (chip is not None): if verbose: print("\n") @@ -74,6 +70,7 @@ class MockObject(object): if verbose: print("After field distortion:\n") print("x = %.2f, y = %.2f\n" % (self.posImg.x, self.posImg.y), flush=True) + x, y = self.posImg.x + 0.5, self.posImg.y + 0.5 self.x_nominal = int(np.floor(x + 0.5)) self.y_nominal = int(np.floor(y + 0.5)) @@ -81,50 +78,22 @@ class MockObject(object): dy = y - self.y_nominal self.offset = galsim.PositionD(dx, dy) + # Deal with chip rotation if chip_wcs is not None: - self.real_wcs = chip_wcs + self.chip_wcs = chip_wcs elif img_header is not None: - self.real_wcs = galsim.FitsWCS(header=img_header) + self.chip_wcs = galsim.FitsWCS(header=img_header) else: - self.real_wcs = None + self.chip_wcs = None - return self.posImg, self.offset, self.localWCS, self.real_wcs, self.fd_shear + return self.posImg, self.offset, self.localWCS, self.chip_wcs, self.fd_shear def getRealPos(self, img, global_x=0., global_y=0., img_real_wcs=None): img_global_pos = galsim.PositionD(global_x, global_y) cel_pos = img.wcs.toWorld(img_global_pos) realPos = img_real_wcs.toImage(cel_pos) - return realPos - def drawObject(self, img, final, flux=None, filt=None, tel=None, exptime=150.): - """ Draw (point like) object on img. - Should be overided for extended source, e.g. galaxy... - Paramter: - img: the "canvas" - final: final (after shear, PSF etc.) GSObject - Return: - img: the image with the GSObject added (or discarded) - isUpdated: is the "canvas" been updated? (a flag for updating statistcs) - """ - isUpdated = True - - # Draw with FFT - # stamp = final.drawImage(wcs=self.localWCS, offset=self.offset) - - # Draw with Photon Shoot - stamp = final.drawImage(wcs=self.localWCS, method='phot', offset=self.offset) - - stamp.setCenter(self.x_nominal, self.y_nominal) - if np.sum(np.isnan(stamp.array)) >= 1: - stamp.setZero() - bounds = stamp.bounds & img.bounds - if bounds.area() == 0: - isUpdated = False - else: - img[bounds] += stamp[bounds] - return img, stamp, isUpdated - def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None): if nphotons_tot == None: @@ -139,29 +108,27 @@ class MockObject(object): self.logger.error(e) return 2, None - nphotons_sum = 0 - photons_list = [] - xmax, ymax = 0, 0 - - # (TEST) Galsim Parameters + # Set Galsim Parameters if self.getMagFilter(filt) <= 15: folding_threshold = 5.e-4 else: folding_threshold = 5.e-3 gsp = galsim.GSParams(folding_threshold=folding_threshold) - + + # Get real image position of object (deal with chip rotation w.r.t its center) self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, - img_real_wcs=self.real_wcs) - + img_real_wcs=self.chip_wcs) x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) y_nominal = int(np.floor(y + 0.5)) dx = x - x_nominal dy = y - y_nominal offset = galsim.PositionD(dx, dy) + # Get real local wcs of object (deal with chip rotation w.r.t its center) + chip_wcs_local = self.chip_wcs.local(self.real_pos) + is_updated = 0 - real_wcs_local = self.real_wcs.local(self.real_pos) - + # Loop over all sub-bandpasses for i in range(len(bandpass_list)): bandpass = bandpass_list[i] try: @@ -170,61 +137,41 @@ class MockObject(object): print(e) if self.logger: self.logger.error(e) - # return False continue - ratio = sub / full - if not (ratio == -1 or (ratio != ratio)): nphotons = ratio * nphotons_tot else: - # return False continue - nphotons_sum += nphotons + + # nphotons_sum += nphotons # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) + + # Get PSF model psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) star = galsim.DeltaFunction(gsparams=gsp) star = star.withFlux(nphotons) star = galsim.Convolve(psf, star) - stamp = star.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) - xmax = max(xmax, stamp.xmax) - ymax = max(ymax, stamp.ymax) - photons = stamp.photons - photons.x += x_nominal - photons.y += y_nominal - photons_list.append(photons) - - stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) - stamp.wcs = real_wcs_local - stamp.setCenter(x_nominal, y_nominal) - bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) + stamp = star.drawImage(wcs=chip_wcs_local, offset=offset) + if np.sum(np.isnan(stamp.array)) > 0: + continue + stamp.setCenter(x_nominal, y_nominal) + bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) + if bounds.area() > 0: + chip.img.setOrigin(0, 0) + chip.img[bounds] += stamp[bounds] + is_updated = 1 + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + del stamp - # # (DEBUG) - # print("stamp bounds: ", stamp.bounds) - # print(bounds) - if bounds.area() > 0: - chip.img.setOrigin(0, 0) - stamp[bounds] = chip.img[bounds] - for i in range(len(photons_list)): - if i == 0: - chip.sensor.accumulate(photons_list[i], stamp) - else: - chip.sensor.accumulate(photons_list[i], stamp, resume=True) - - chip.img[bounds] = stamp[bounds] - - chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) - else: - # Return code 0: object photons missed this detector + if is_updated == 0: + # Return code 0: object has missed this detector print("obj %s missed"%(self.id)) if self.logger: self.logger.info("obj %s missed"%(self.id)) return 0, pos_shear - - del photons_list - del stamp return 1, pos_shear # Return code 1: draw sucesss def addSLStoChipImage(self, sdp=None, chip=None, xOrderSigPlus=None, local_wcs=None): @@ -299,7 +246,7 @@ class MockObject(object): names=('WAVELENGTH', 'FLUX')) self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, - img_real_wcs=self.real_wcs) + img_real_wcs=self.chip_wcs) x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) @@ -308,7 +255,7 @@ class MockObject(object): dy = y - y_nominal offset = galsim.PositionD(dx, dy) - real_wcs_local = self.real_wcs.local(self.real_pos) + chip_wcs_local = self.chip_wcs.local(self.real_pos) flat_cube = chip.flat_cube @@ -323,7 +270,7 @@ class MockObject(object): star = star.withFlux(tel.pupil_area * exptime) star = galsim.Convolve(psf, star) - starImg = star.drawImage(nx=100, ny=100, wcs=real_wcs_local, offset=offset) + starImg = star.drawImage(nx=100, ny=100, wcs=chip_wcs_local, offset=offset) origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_nominal - (starImg.center.x - starImg.xmin)] @@ -350,7 +297,7 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) @@ -367,7 +314,7 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) del sdp_p1 del sdp_p2 @@ -379,7 +326,7 @@ class MockObject(object): conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) del sdp elif grating_split_pos_chip >= gal_end[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, @@ -389,7 +336,7 @@ class MockObject(object): conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) del sdp del psf return 1, pos_shear diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index a70f997..f3001ac 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -186,7 +186,7 @@ class Observation(object): extName='SCI', timestamp = pointing.timestamp, exptime = pointing.exp_time, - readoutTime = 40.) + readoutTime = chip.readout_time) chip_wcs = galsim.FitsWCS(header=h_ext) @@ -382,7 +382,7 @@ class Observation(object): extName='SCI', timestamp=pointing.timestamp, exptime=pointing.exp_time, - readoutTime=40.) + readoutTime=chip.readout_time) chip.img = galsim.Image(chip.img.array, dtype=np.uint16) hdu1 = fits.PrimaryHDU(header=h_prim) diff --git a/config/config_50sqdeg.yaml b/config/config_50sqdeg.yaml index c0f4e7f..955d878 100644 --- a/config/config_50sqdeg.yaml +++ b/config/config_50sqdeg.yaml @@ -9,9 +9,9 @@ # Base diretories and naming setup # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent -work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/new_sim/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "rotate_0" +run_name: "test_new_sim" # Whether to use MPI run_option: @@ -45,10 +45,10 @@ catalog_options: # AGN_SED_WAVE: "wave_ross13.npy" # Only simulate stars? - star_only: NO + star_only: YES # Only simulate galaxies? - galaxy_only: YES + galaxy_only: NO # rotate galaxy ellipticity rotateEll: 0. # [degree] @@ -171,7 +171,7 @@ ins_effects: non_linear: ON # Whether to add non-linearity cosmic_ray: ON # Whether to add cosmic-ray cray_differ: ON # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: ON # Whether to simulate CTE trails + cte_trail: OFF # Whether to simulate CTE trails saturbloom: ON # Whether to simulate Saturation & Blooming add_badcolumns: ON # Whether to add bad columns add_hotpixels: ON # Whether to add hot pixels diff --git a/libmoduleDfBF/diffusion_X1.c b/libmoduleDfBF/diffusion_X1.c new file mode 100644 index 0000000..199b48c --- /dev/null +++ b/libmoduleDfBF/diffusion_X1.c @@ -0,0 +1,120 @@ +#include +#include +#include +#include +#include "nrutil.h" + +#define ISSETBITFLAG(x,b) ((x) & (1 << (b))) +#define ADD_DIFFUSION 1 +#define ADD_BF_FILTER 2 + +float linearInterp(float xp, float x0, float y0, float x1, float y1); + +void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bit_flag) +{ + int nx, ny, i,j,k,ks; + int it,jt,itt,jtt; + int diffuidx[26][2],diffuN,ilow,ih,im,dim[3]; + float diffua[5][5],cdiffu[26],**bfa; + double mvar,mcov,tmp,ma,mb,mc; + char fname[100]; + + nx = ngx_ima; //input-image size + ny = ngy_ima; + + //0. init. original image with an input array (arr_ima) + //1. Adding diffusion effect. + if(ISSETBITFLAG(bit_flag, ADD_DIFFUSION)) + { + printf("adding diffusion.....\n"); + printf("ERR: no diffusion filter ..."); + exit(0); + } + + + //2. Adding BF effect + if(ISSETBITFLAG(bit_flag, ADD_BF_FILTER)) + { + printf("Adding BF effect...\n"); + //setup BF correlation fliter + float neX; + float neP1 = 50000; + float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302}; + float neP2 = 10000; + float bfaP2[9]={0.9945288, 0.0003041936, 0.0007539311, 0.0002424675, 0.0001226098, 0.00009308617, 0.00008027447, 0.00006309676, 0.00006400052}; + + bfa=matrix(-2,2,-2,2); + + // smooth with the BF filter + for(i=0;i= 10000) + { + bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0; + bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575; + bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652; + bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335; + bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]); + bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118; + bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]); + bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083; + bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043; + } + else + { + neX=10000; + bfa[0][0]=0; + bfa[0][1]=bfa[0][-1]=bfaP2[1]; + bfa[-1][0]=bfa[1][0]=bfaP2[2]; + bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=bfaP2[3]; + bfa[0][-2]=bfa[0][2]=bfaP2[4]; + bfa[-2][0]=bfa[2][0]=bfaP2[5]; + bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=bfaP2[6]; + bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=bfaP2[7]; + bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=bfaP2[8]; + } + tmp = 0; + for(it=-2;it<=2;it++) + for(jt=-2;jt<=2;jt++) + { + bfa[it][jt] = bfa[it][jt]/neX*arr_ima[j+i*ny]; + tmp += bfa[it][jt]; + } + bfa[0][0]=1.-tmp; + + // assign electrons according to the BF filter bfat + for(it=-2;it<=2;it++) + { + for(jt=-2;jt<=2;jt++) + { + itt=i+it; + jtt=j+jt; + if(itt>=0 && jtt>=0 && itt