diff --git a/ObservationSim/Config/Config.py b/ObservationSim/Config/Config.py index 74837df4fcc5e42ff506a39a9dae2ab572c4d53b..3d6724b139c2fc2f6cad907185d4947e04d34329 100755 --- a/ObservationSim/Config/Config.py +++ b/ObservationSim/Config/Config.py @@ -20,6 +20,7 @@ def config_dir(config, work_dir=None, data_dir=None): # PSF data directory if config["psf_setting"]["psf_model"] == "Interp": path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"]) + path_dict["psf_sls_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_sls_dir"]) return path_dict diff --git a/ObservationSim/Instrument/Chip/ChipUtils.py b/ObservationSim/Instrument/Chip/ChipUtils.py index d8e9a0aaf4cfb895f1e173599a2015bd95688a93..99b5ddce4ca392d25f054151f87a49b1840446b4 100644 --- a/ObservationSim/Instrument/Chip/ChipUtils.py +++ b/ObservationSim/Instrument/Chip/ChipUtils.py @@ -21,6 +21,22 @@ def log_info(msg, logger=None): else: print(msg, flush=True) +def getChipSLSGratingID(chipID): + gratingID = ['',''] + if chipID == 1: gratingID = ['GI2', 'GI1'] + if chipID == 2: gratingID = ['GV4', 'GV3'] + if chipID == 3: gratingID = ['GU2', 'GU1'] + if chipID == 4: gratingID = ['GU4', 'GU3'] + if chipID == 5: gratingID = ['GV2', 'GV1'] + if chipID == 10: gratingID = ['GI4', 'GI3'] + if chipID == 21: gratingID = ['GI6', 'GI5'] + if chipID == 26: gratingID = ['GV8', 'GV7'] + if chipID == 27: gratingID = ['GU6', 'GU5'] + if chipID == 28: gratingID = ['GU8', 'GU7'] + if chipID == 29: gratingID = ['GV6', 'GV5'] + if chipID == 30: gratingID = ['GI8', 'GI7'] + return gratingID + def getChipSLSConf(chipID): confFile = '' if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf'] diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 81ce517cb7bc4fbcefea1905dac16afaf63e2322..e7f43216a7a71a451fbd1815573b1f33cc7c51f3 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -248,10 +248,27 @@ class Galaxy(MockObject): xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} grating_split_pos_chip = 0 + grating_split_pos + + branges = np.zeros([len(bandpass_list), 2]) + + # print(hasattr(psf_model, 'bandranges')) + + if hasattr(psf_model, 'bandranges'): + if psf_model.bandranges is None: + return 2, None + if len(psf_model.bandranges) != len(bandpass_list): + return 2, None + branges = psf_model.bandranges + else: + for i in range(len(bandpass_list)): + branges[i, 0] = bandpass_list[i].blue_limit * 10 + branges[i, 1] = bandpass_list[i].red_limit * 10 + for i in range(len(bandpass_list)): - bandpass = bandpass_list[i] + # bandpass = bandpass_list[i] + brange = branges[i] - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) + # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) 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) @@ -272,14 +289,14 @@ class Galaxy(MockObject): g2 += fd_shear.g2 gal_shear = galsim.Shear(g1=g1, g2=g2) gal = gal.shear(gal_shear) - gal = galsim.Convolve(psf, gal) + # gal = galsim.Convolve(psf, gal) - if not big_galaxy: # Not apply PSF for very big galaxy - gal = galsim.Convolve(psf, gal) - # if fd_shear is not None: - # gal = gal.shear(fd_shear) + # if not big_galaxy: # Not apply PSF for very big galaxy + # gal = galsim.Convolve(psf, gal) + # # if fd_shear is not None: + # # gal = gal.shear(fd_shear) - starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset) + starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space') origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_nominal - (starImg.center.x - starImg.xmin)] @@ -301,12 +318,16 @@ class Galaxy(MockObject): sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local) subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) @@ -318,12 +339,16 @@ class Galaxy(MockObject): sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local) del sdp_p1 del sdp_p2 @@ -331,25 +356,33 @@ class Galaxy(MockObject): sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local) del sdp elif grating_split_pos_chip>=gal_end[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) - del psf + # del psf return 1, pos_shear def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.): diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 654cb720188cd71a1b98e541cac26e866dec2160..0dd8fb0543eb5efef5373b520c3059c5fb628e2c 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -4,7 +4,7 @@ import astropy.constants as cons from astropy import wcs from astropy.table import Table -from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders +from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders, convolveImg from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, \ getABMAG from ObservationSim.MockObject.SpecDisperser import SpecDisperser @@ -222,9 +222,69 @@ class MockObject(object): del stamp del spec_orders + def addSLStoChipImageWithPSF(self, sdp=None, chip=None, pos_img_local = [1,1], psf_model=None, bandNo = 1, grating_split_pos=3685, local_wcs=None, pos_img=None): + spec_orders = sdp.compute_spec_orders() + for k, v in spec_orders.items(): + img_s = v[0] + # print(bandNo,k) + try: + psf, pos_shear = psf_model.get_PSF(chip, pos_img_local = pos_img_local, bandNo = bandNo, galsimGSObject=True, g_order = k, grating_split_pos=grating_split_pos) + except: + psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) + + psf_img = psf.drawImage(nx=100, ny=100, wcs = local_wcs) + + psf_img_m = psf_img.array + + ######################################################### + # DEBUG + ######################################################### + # ids_p = psf_img_m < 0 + # psf_img_m[ids_p] = 0 + + # from astropy.io import fits + # fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m) + + # print("DEBUG: orig_off is", orig_off) + nan_ids = np.isnan(img_s) + if img_s[nan_ids].shape[0] > 0: + img_s[nan_ids] = 0 + print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0]) + ######################################################### + img_s, orig_off = convolveImg(img_s, psf_img_m) + origin_order_x = v[1] - orig_off[0] + origin_order_y = v[2] - orig_off[1] + + + specImg = galsim.ImageF(img_s) + # photons = galsim.PhotonArray.makeFromImage(specImg) + # photons.x += origin_order_x + # photons.y += origin_order_y + + # xlen_imf = int(specImg.xmax - specImg.xmin + 1) + # ylen_imf = int(specImg.ymax - specImg.ymin + 1) + # stamp = galsim.ImageF(xlen_imf, ylen_imf) + # stamp.wcs = local_wcs + # stamp.setOrigin(origin_order_x, origin_order_y) + + specImg.wcs = local_wcs + specImg.setOrigin(origin_order_x, origin_order_y) + + bounds = specImg.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) + if bounds.area() == 0: + continue + chip.img.setOrigin(0, 0) + chip.img[bounds] = chip.img[bounds] + specImg[bounds] + # stamp[bounds] = chip.img[bounds] + # # chip.sensor.accumulate(photons, stamp) + # chip.img[bounds] = stamp[bounds] + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + # del stamp + del spec_orders + return pos_shear + def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None): - if normFilter is not None: norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed, @@ -262,22 +322,38 @@ class MockObject(object): xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062, 'C': 4.035447379743442, 'D': 5.5684364343742825, 'E': 16.260021029735388} grating_split_pos_chip = 0 + grating_split_pos + + branges = np.zeros([len(bandpass_list),2]) + + if hasattr(psf_model,'bandranges'): + if psf_model.bandranges is None: + return 2, None + if len(psf_model.bandranges) != len(bandpass_list): + return 2, None + branges = psf_model.bandranges + else: + for i in range(len(bandpass_list)): + branges[i, 0] = bandpass_list[i].blue_limit * 10 + branges[i, 1] = bandpass_list[i].red_limit * 10 + for i in range(len(bandpass_list)): - bandpass = bandpass_list[i] - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, - folding_threshold=folding_threshold) + # bandpass = bandpass_list[i] + brange = branges[i] + + # 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(tel.pupil_area * exptime) - star = galsim.Convolve(psf, star) + psf_tmp = galsim.Gaussian(sigma=0.002) + star = galsim.Convolve(psf_tmp, star) - starImg = star.drawImage(nx=100, ny=100, wcs=chip_wcs_local, offset=offset) + starImg = star.drawImage(nx=60, ny=60, wcs=chip_wcs_local, offset=offset) origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_nominal - (starImg.center.x - starImg.xmin)] starImg.setOrigin(0,0) gal_origin = [origin_star[0], origin_star[1]] gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1] - if gal_origin[1] < grating_split_pos_chip < gal_end[1]: subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) ## part img disperse @@ -292,12 +368,15 @@ class MockObject(object): sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local = [xcenter_p1,ycenter_p1], + psf_model=psf_model, bandNo = i+1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) @@ -309,12 +388,15 @@ class MockObject(object): sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) del sdp_p1 del sdp_p2 @@ -322,23 +404,29 @@ class MockObject(object): sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) del sdp elif grating_split_pos_chip >= gal_end[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, - band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) del sdp - del psf + # del psf return 1, pos_shear def SNRestimate(self, img_obj, flux, noise_level=0.0, seed=31415): diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py index 3f3f52b01a1d9b061099491d5e49f677979411ee..7c8b51efa9f39b07e485f00506a53eaff163f118 100644 --- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +++ b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py @@ -155,7 +155,7 @@ class SpecDisperser(object): sensitivity_beam = ysens len_spec_x = len(dx) - len_spec_y = int(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0]) + 1) + len_spec_y = int(abs(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0])) + 1) beam_sh = (self.img_sh[0] + len_spec_y, self.img_sh[1] + len_spec_x) modelf = zeros(product(beam_sh), dtype=float) diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py index 6f51570bd9426fd40c94b019f17914cde6f65150..02aee7700f411d53bf05850bfee514d3f4898f45 100755 --- a/ObservationSim/MockObject/_util.py +++ b/ObservationSim/MockObject/_util.py @@ -571,4 +571,14 @@ def convolveGaussXorders(img=None, sigma = 1): convImg = signal.fftconvolve(img, psf, mode='full', axes=None) return convImg, offset +def convolveImg(img=None, psf = None): + from astropy.modeling.models import Gaussian2D + from scipy import signal + + convImg = signal.fftconvolve(img, psf, mode='full', axes=None) + offset_x = int(psf.shape[1]/2. + 0.5) - 1 + offset_y = int(psf.shape[0]/2. + 0.5) - 1 + offset = [offset_x,offset_y] + return convImg, offset + diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index f3001ac1b4c117dc883ae227a11088892c77d674..d1e10f76a5eeaff0697c56155c9873d471a51256 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -15,7 +15,7 @@ from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensio from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip from ObservationSim.Instrument.Chip import Effects from ObservationSim.Straylight import calculateSkyMap_split_g -from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp +from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS from ObservationSim._util import get_shear_field, makeSubDir_PointingList from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position @@ -63,7 +63,10 @@ class Observation(object): if self.config["psf_setting"]["psf_model"] == "Gauss": psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) elif self.config["psf_setting"]["psf_model"] == "Interp": - psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"]) + if chip.survey_type == "spectroscopic": + psf_model = PSFInterpSLS(chip, filt,PSF_data_prefix=self.path_dict["psf_sls_dir"]) + else: + psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"]) else: chip_output.Log_error("unrecognized PSF model type!!", flush=True) @@ -198,6 +201,10 @@ class Observation(object): obj = self.cat.objs[j] + # (DEBUG) + # if obj.getMagFilter(filt)>20: + # continue + # load and convert SED; also caculate object's magnitude in all CSST bands try: sed_data = self.cat.load_sed(obj) diff --git a/ObservationSim/PSF/PSFInterpSLS.py b/ObservationSim/PSF/PSFInterpSLS.py new file mode 100644 index 0000000000000000000000000000000000000000..7dc1624c2ae4305cac78650044e46c95e6de57a2 --- /dev/null +++ b/ObservationSim/PSF/PSFInterpSLS.py @@ -0,0 +1,570 @@ +''' +PSF interpolation for CSST-Sim + +NOTE: [iccd, iwave, ipsf] are counted from 1 to n, but [tccd, twave, tpsf] are counted from 0 to n-1 +''' + +import sys +import time +import copy +import numpy as np +import scipy.spatial as spatial +import galsim +import h5py + +from ObservationSim.PSF.PSFModel import PSFModel +from ObservationSim.Instrument.Chip import ChipUtils as chip_utils +import os +from astropy.io import fits + +from astropy.modeling.models import Gaussian2D +from scipy import signal + + +LOG_DEBUG = False #***# +NPSF = 900 #***# 30*30 +PixSizeInMicrons = 5. #***# in microns + + +###find neighbors-KDtree### +# def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True): +# """ +# find nearest neighbors by 2D-KDTree +# +# Parameters: +# tx, ty (float, float): a given position +# px, py (numpy.array, numpy.array): position data for tree +# dr (float-optional): distance +# dn (int-optional): nearest-N +# OnlyDistance (bool-optional): only use distance to find neighbors. Default: True +# +# Returns: +# dataq (numpy.array): index +# """ +# datax = px +# datay = py +# tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel()))) +# +# dataq=[] +# rr = dr +# if OnlyDistance == True: +# dataq = tree.query_ball_point([tx, ty], rr) +# if OnlyDistance == False: +# while len(dataq) < dn: +# dataq = tree.query_ball_point([tx, ty], rr) +# rr += dr +# dd = np.hypot(datax[dataq]-tx, datay[dataq]-ty) +# ddSortindx = np.argsort(dd) +# dataq = np.array(dataq)[ddSortindx[0:dn]] +# return dataq +# +# ###find neighbors-hoclist### +# def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy): +# if np.max(partx) > nhocx*dhocx: +# print('ERROR') +# sys.exit() +# if np.max(party) > nhocy*dhocy: +# print('ERROR') +# sys.exit() +# +# npart = partx.size +# hoclist= np.zeros(npart, dtype=np.int32)-1 +# hoc = np.zeros([nhocy, nhocx], dtype=np.int32)-1 +# for ipart in range(npart): +# ix = int(partx[ipart]/dhocx) +# iy = int(party[ipart]/dhocy) +# hoclist[ipart] = hoc[iy, ix] +# hoc[iy, ix] = ipart +# return hoc, hoclist +# +# def hocFind(px, py, dhocx, dhocy, hoc, hoclist): +# ix = int(px/dhocx) +# iy = int(py/dhocy) +# +# neigh=[] +# it = hoc[iy, ix] +# while it != -1: +# neigh.append(it) +# it = hoclist[it] +# return neigh +# +# def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None): +# nhocy = nhocx = 20 +# +# pxMin = np.min(px) +# pxMax = np.max(px) +# pyMin = np.min(py) +# pyMax = np.max(py) +# +# dhocx = (pxMax - pxMin)/(nhocx-1) +# dhocy = (pyMax - pyMin)/(nhocy-1) +# partx = px - pxMin +dhocx/2 +# party = py - pyMin +dhocy/2 +# +# if hoc is None: +# hoc, hoclist = hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy) +# return hoc, hoclist +# +# if hoc is not None: +# tx = tx - pxMin +dhocx/2 +# ty = ty - pyMin +dhocy/2 +# itx = int(tx/dhocx) +# ity = int(ty/dhocy) +# +# ps = [-1, 0, 1] +# neigh=[] +# for ii in range(3): +# for jj in range(3): +# ix = itx + ps[ii] +# iy = ity + ps[jj] +# if ix < 0: +# continue +# if iy < 0: +# continue +# if ix > nhocx-1: +# continue +# if iy > nhocy-1: +# continue +# +# #neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist) +# it = hoc[iy, ix] +# while it != -1: +# neigh.append(it) +# it = hoclist[it] +# #neigh.append(neightt) +# #ll = [i for k in neigh for i in k] +# if dn != -1: +# ptx = np.array(partx[neigh]) +# pty = np.array(party[neigh]) +# dd = np.hypot(ptx-tx, pty-ty) +# idx = np.argsort(dd) +# neigh= np.array(neigh)[idx[0:dn]] +# return neigh +# +# +# ###PSF-IDW### +# def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False): +# """ +# psf interpolation by IDW +# +# Parameters: +# px, py (float, float): position of the target +# PSFMat (numpy.array): image +# cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers +# IDWindex (int-optional): the power index of IDW +# OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation +# +# Returns: +# psfMaker (numpy.array) +# """ +# +# minimum_psf_weight = 1e-8 +# ref_col = px +# ref_row = py +# +# ngy, ngx = PSFMat[0, :, :].shape +# npsf = PSFMat[:, :, :].shape[0] +# psfWeight = np.zeros([npsf]) +# +# if OnlyNeighbors == True: +# if hoc is None: +# neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False) +# if hoc is not None: +# neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist) +# +# neighFlag = np.zeros(npsf) +# neighFlag[neigh] = 1 +# +# for ipsf in range(npsf): +# if OnlyNeighbors == True: +# if neighFlag[ipsf] != 1: +# continue +# +# dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2) +# if IDWindex == 1: +# psfWeight[ipsf] = dist +# if IDWindex == 2: +# psfWeight[ipsf] = dist**2 +# if IDWindex == 3: +# psfWeight[ipsf] = dist**3 +# if IDWindex == 4: +# psfWeight[ipsf] = dist**4 +# psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight) +# psfWeight[ipsf] = 1./psfWeight[ipsf] +# psfWeight /= np.sum(psfWeight) +# +# psfMaker = np.zeros([ngy, ngx], dtype=np.float32) +# for ipsf in range(npsf): +# if OnlyNeighbors == True: +# if neighFlag[ipsf] != 1: +# continue +# +# iPSFMat = PSFMat[ipsf, :, :].copy() +# ipsfWeight = psfWeight[ipsf] +# +# psfMaker += iPSFMat * ipsfWeight +# psfMaker /= np.nansum(psfMaker) +# +# return psfMaker + + + +###define PSFInterp### +class PSFInterpSLS(PSFModel): + def __init__(self, chip, filt,PSF_data_prefix="", sigSpin=0, psfRa=0.15, pix_size = 0.005): + if LOG_DEBUG: + print('===================================================') + print('DEBUG: psf module for csstSim ' \ + +time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True) + print('===================================================') + + self.sigSpin = sigSpin + self.sigGauss = psfRa + self.grating_ids = chip_utils.getChipSLSGratingID(chip.chipID) + _,self.grating_type = chip.getChipFilter(chipID=chip.chipID) + self.data_folder = PSF_data_prefix + self.getPSFDataFromFile(filt) + self.pixsize = pix_size # um + + def getPSFDataFromFile(self, filt): + gratingInwavelist = {'GU':0,'GV':1,'GI':2} + grating_orders = ['0','1'] + waveListFn = self.data_folder + '/wavelist.dat' + wavelists = np.loadtxt(waveListFn) + self.waveList = wavelists[:,gratingInwavelist[self.grating_type]] + bandranges = np.zeros([4,2]) + midBand = (self.waveList[0:3] + self.waveList[1:4])/2.*10000. + bandranges[0,0] = filt.blue_limit + bandranges[1:4,0] = midBand + bandranges[0:3, 1] = midBand + bandranges[3,1] = filt.red_limit + + self.bandranges = bandranges + + self.grating1_data = {} + g_folder = self.data_folder + '/' + self.grating_ids[0] + '/' + for g_order in grating_orders: + g_folder_order = g_folder + 'PSF_Order_' + g_order + '/' + grating_order_data = {} + for bandi in [1,2,3,4]: + subBand_data = {} + subBand_data['bandrange'] = bandranges[bandi-1] + final_folder = g_folder_order + str(bandi) + '/' + print(final_folder) + pca_fs = os.listdir(final_folder) + for fname in pca_fs: + if ('_PCs.fits' in fname) and (fname[0] != '.'): + fname_ = final_folder + fname + hdu = fits.open(fname_) + subBand_data['band_data'] = hdu + grating_order_data['band'+str(bandi)] = subBand_data + self.grating1_data['order'+g_order] = grating_order_data + + self.grating2_data = {} + g_folder = self.data_folder + '/' + self.grating_ids[1] + '/' + for g_order in grating_orders: + g_folder_order = g_folder + 'PSF_Order_' + g_order + '/' + grating_order_data = {} + for bandi in [1, 2, 3, 4]: + subBand_data = {} + subBand_data['bandrange'] = bandranges[bandi - 1] + final_folder = g_folder_order + str(bandi) + '/' + print(final_folder) + pca_fs = os.listdir(final_folder) + for fname in pca_fs: + if ('_PCs.fits' in fname) and (fname[0] != '.'): + fname_ = final_folder + fname + hdu = fits.open(fname_) + subBand_data['band_data'] = hdu + grating_order_data['band' + str(bandi)] = subBand_data + self.grating2_data['order' + g_order] = grating_order_data + + # + # + # + # def _getPSFwave(self, iccd, PSF_data_file, PSF_data_prefix): + # # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r') + # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r') + # nwave = len(fq.keys()) + # fq.close() + # return nwave + # + # + # def _loadPSF(self, iccd, PSF_data_file, PSF_data_prefix): + # psfSet = [] + # # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r') + # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r') + # for ii in range(self.nwave): + # iwave = ii+1 + # psfWave = [] + # + # fq_iwave = fq['w_{:}'.format(iwave)] + # for jj in range(self.npsf): + # ipsf = jj+1 + # psfInfo = {} + # psfInfo['wavelength']= fq_iwave['wavelength'][()] + # + # fq_iwave_ipsf = fq_iwave['psf_{:}'.format(ipsf)] + # psfInfo['pixsize'] = PixSizeInMicrons + # psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()] + # psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()] + # psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()] + # psfInfo['image_y'] = fq_iwave_ipsf['image_y'][()] + # psfInfo['centroid_x']= fq_iwave_ipsf['cx'][()] + # psfInfo['centroid_y']= fq_iwave_ipsf['cy'][()] + # psfInfo['psfMat'] = fq_iwave_ipsf['psfMat'][()] + # + # psfWave.append(psfInfo) + # psfSet.append(psfWave) + # fq.close() + # + # if LOG_DEBUG: + # print('psfSet has been loaded:', flush=True) + # print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True) + # return psfSet + # + # + # def _findWave(self, bandpass): + # if isinstance(bandpass,int): + # twave = bandpass + # return twave + # + # for twave in range(self.nwave): + # bandwave = self.PSF_data[twave][0]['wavelength'] + # if bandpass.blue_limit < bandwave and bandwave < bandpass.red_limit: + # return twave + # return -1 + # + # + + def convolveWithGauss(self, img=None, sigma=1): + + offset = int(np.ceil(sigma * 3)) + g_size = 2 * offset + 1 + m_cen = int(g_size / 2) + print('-----',g_size) + g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma) + yp, xp = np.mgrid[0:g_size, 0:g_size] + g_PSF = g_PSF_(xp, yp) + psf = g_PSF / g_PSF.sum() + convImg = signal.fftconvolve(img, psf, mode='full', axes=None) + convImg = convImg/np.sum(convImg) + return convImg + + + def get_PSF(self, chip, pos_img_local = [1000,1000], bandNo = 1, galsimGSObject=True, folding_threshold=5.e-3, g_order = 'A', grating_split_pos=3685): + """ + Get the PSF at a given image position + + Parameters: + chip: A 'Chip' object representing the chip we want to extract PSF from. + pos_img: A 'galsim.Position' object representing the image position. + bandpass: A 'galsim.Bandpass' object representing the wavelength range. + pixSize: The pixels size of psf matrix + findNeighMode: 'treeFind' or 'hoclistFind' + Returns: + PSF: A 'galsim.GSObject'. + """ + order_IDs = {'A': '1', 'B': '0' ,'C': '0', 'D': '0', 'E': '0'} + contam_order_sigma = {'C':0.28032344707964174,'D':0.39900182912061344,'E':1.1988309797685412} #arcsec + x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. + y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. + # print(pos_img.x - x_start) + pos_img_x = pos_img_local[0] + x_start + pos_img_y = pos_img_local[1] + y_start + pos_img = galsim.PositionD(pos_img_x, pos_img_y) + if pos_img_local[0] < grating_split_pos: + psf_data = self.grating1_data + else: + psf_data = self.grating2_data + + grating_order = order_IDs[g_order] + # if grating_order in ['-2','-1','2']: + # grating_order = '1' + + + # if grating_order in ['0', '1']: + psf_order = psf_data['order'+grating_order] + psf_order_b = psf_order['band'+str(bandNo)] + psf_b_dat = psf_order_b['band_data'] + pos_p = psf_b_dat[1].data + pc_coeff = psf_b_dat[2].data + pcs = psf_b_dat[0].data + # print(max(pos_p[:,0]), min(pos_p[:,0]),max(pos_p[:,1]), min(pos_p[:,1])) + # print(chip.x_cen, chip.y_cen) + # print(pos_p) + px = pos_img.x*chip.pix_size + py = pos_img.y*chip.pix_size + + dist2=(pos_p[:,1] - px)*(pos_p[:,1] - px) + (pos_p[:,0] - py)*(pos_p[:,0] - py) + temp_sort_dist = np.zeros([dist2.shape[0],2]) + temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0],1) + temp_sort_dist[:, 1] = dist2 + # print(temp_sort_dist) + dits2_sortlist = sorted(temp_sort_dist, key=lambda x:x[1]) + # print(dits2_sortlist) + nearest4p = np.zeros([4,2]) + pc_coeff_4p = np.zeros([pc_coeff.data.shape[0],4]) + + for i in np.arange(4): + smaller_ids = int(dits2_sortlist[i][0]) + nearest4p[i, 0] = pos_p[smaller_ids, 1] + nearest4p[i, 1] = pos_p[smaller_ids, 0] + pc_coeff_4p[:,i] = pc_coeff[:,smaller_ids] + idw_dist = 1/(np.sqrt((px-nearest4p[:,0]) * (px-nearest4p[:,0]) + (py-nearest4p[:,1]) * (py-nearest4p[:,1]))) + + coeff_int = np.zeros(pc_coeff.data.shape[0]) + for i in np.arange(4): + coeff_int = coeff_int + pc_coeff_4p[:,i]*idw_dist[i] + coeff_int = coeff_int / np.sum(coeff_int) + + npc = 10 + m_size = int(pcs.shape[0]**0.5) + PSF_int = np.dot(pcs[:,0:npc],coeff_int[0:npc]).reshape(m_size,m_size) + + # PSF_int = PSF_int/np.sum(PSF_int) + PSF_int_trans = np.flipud(np.fliplr(PSF_int)) + PSF_int_trans = np.fliplr(PSF_int_trans.T) + # PSF_int_trans = np.abs(PSF_int_trans) + # ids_szero = PSF_int_trans<0 + # PSF_int_trans[ids_szero] = 0 + # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape) + PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans) + + # from astropy.io import fits + # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans) + + + # if g_order in ['C','D','E']: + # g_simgma = contam_order_sigma[g_order]/pixel_size_arc + # PSF_int_trans = self.convolveWithGauss(PSF_int_trans,g_simgma) + # n_m_size = int(m_size/2) + # + # n_PSF_int = np.zeros([n_m_size, n_m_size]) + # + # for i in np.arange(n_m_size): + # for j in np.arange(n_m_size): + # n_PSF_int[i,j] = np.sum(PSF_int[2*i:2*i+2, 2*j:2*j+2]) + # + # n_PSF_int = n_PSF_int/np.sum(n_PSF_int) + + # chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) + # chip.img.wcs = galsim.wcs.AffineTransform + if galsimGSObject: + + # imPSFt = np.zeros([257,257]) + # imPSFt[0:256, 0:256] = imPSF + # # imPSFt[120:130, 0:256] = 1. + + pixel_size_arc = np.rad2deg(self.pixsize * 1e-3 / 28) * 3600 + img = galsim.ImageF(PSF_int_trans, scale=pixel_size_arc) + gsp = galsim.GSParams(folding_threshold=folding_threshold) + ############TEST: START + # Use sheared PSF to test the PSF orientation + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + ############TEST: END + self.psf = galsim.InterpolatedImage(img, gsparams=gsp) + # if g_order in ['C','D','E']: + # add_psf = galsim.Gaussian(sigma=contam_order_sigma[g_order], flux=1.0) + # self.psf = galsim.Convolve(self.psf, add_psf) + wcs = chip.img.wcs.local(pos_img) + scale = galsim.PixelScale(0.074) + self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img)) + + # return self.PSFspin(x=px/0.01, y=py/0.01) + return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) + + return PSF_int_trans, PSF_int + + + # pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize + # + # # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID' + # twave = self._findWave(bandpass) + # if twave == -1: + # print("!!!PSF bandpass does not match.") + # exit() + # PSFMat = self.psfMat[twave] + # cen_col= self.cen_col[twave] + # cen_row= self.cen_row[twave] + # + # px = (pos_img.x - chip.cen_pix_x)*0.01 + # py = (pos_img.y - chip.cen_pix_y)*0.01 + # if findNeighMode == 'treeFind': + # imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, PSFCentroidWgt=True) + # if findNeighMode == 'hoclistFind': + # assert(self.hoc != 0), 'hoclist should be built correctly!' + # imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True) + # + # ############TEST: START + # TestGaussian = False + # if TestGaussian: + # gsx = galsim.Gaussian(sigma=0.04) + # #pointing_pa = -23.433333 + # imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array + # ############TEST: END + # + # if galsimGSObject: + # imPSFt = np.zeros([257,257]) + # imPSFt[0:256, 0:256] = imPSF + # # imPSFt[120:130, 0:256] = 1. + # + # img = galsim.ImageF(imPSFt, scale=pixSize) + # gsp = galsim.GSParams(folding_threshold=folding_threshold) + # ############TEST: START + # # Use sheared PSF to test the PSF orientation + # # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + # ############TEST: END + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp) + # wcs = chip.img.wcs.local(pos_img) + # scale = galsim.PixelScale(0.074) + # self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img)) + # + # # return self.PSFspin(x=px/0.01, y=py/0.01) + # return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) + # return imPSF + # + # def PSFspin(self, x, y): + # """ + # The PSF profile at a given image position relative to the axis center + # + # Parameters: + # theta : spin angles in a given exposure in unit of [arcsecond] + # dx, dy: relative position to the axis center in unit of [pixels] + # + # Return: + # Spinned PSF: g1, g2 and axis ratio 'a/b' + # """ + # a2Rad = np.pi/(60.0*60.0*180.0) + # + # ff = self.sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels] + # rc = np.sqrt(x*x + y*y) + # cpix = rc*(self.sigSpin*a2Rad) + # + # beta = (np.arctan2(y,x) + np.pi/2) + # ell = cpix**2/(2.0*ff**2+cpix**2) + # qr = np.sqrt((1.0+ell)/(1.0-ell)) + # PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians) + # return self.psf.shear(PSFshear), PSFshear + +from ObservationSim.Instrument import Filter, FilterParam, Chip +import yaml +if __name__ == '__main__': + configfn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/config/config_C6_dev.yaml' + with open(configfn, "r") as stream: + try: + config = yaml.safe_load(stream) + for key, value in config.items(): + print (key + " : " + str(value)) + except yaml.YAMLError as exc: + print(exc) + chip = Chip(chipID=1,config=config) + filter_id, filter_type = chip.getChipFilter() + filt = Filter(filter_id=filter_id, + filter_type=filter_type, + filter_param=FilterParam()) + + psf_i = PSFInterpSLS(chip, filt,PSF_data_prefix="/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/") + pos_img = galsim.PositionD(x=25155, y=-22060) + psf_im = psf_i.get_PSF(chip, pos_img = pos_img, g_order = '1') + diff --git a/ObservationSim/PSF/__init__.py b/ObservationSim/PSF/__init__.py index d0ea11850ebb2728337a2b2ce76d714b4591fb5a..80bf7d3685fc1aba1b88aa8a4535163ca39d55fa 100755 --- a/ObservationSim/PSF/__init__.py +++ b/ObservationSim/PSF/__init__.py @@ -2,4 +2,5 @@ from .PSFModel import PSFModel from .PSFGauss import PSFGauss # from .PSFInterp.PSFInterp import PSFInterp from .PSFInterp import PSFInterp +from .PSFInterpSLS import PSFInterpSLS from .FieldDistortion import FieldDistortion \ No newline at end of file diff --git a/config/config_C6_dev.yaml b/config/config_C6_dev.yaml new file mode 100644 index 0000000000000000000000000000000000000000..c939acc95712125fe1639c01e8a1527758606454 --- /dev/null +++ b/config/config_C6_dev.yaml @@ -0,0 +1,221 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# 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: "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/" +data_dir: "/Volumes/EAGET/C6_data/inputData/" +run_name: "profile_C6" + +# Whether to use MPI +run_option: + use_mpi: NO + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "Catalog_C6_20221212" + star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + galaxy_cat: "cat2CSSTSim_bundle/" + AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "quickspeclib_ross13.fits" + AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: NO + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + seed_Av: 121212 # Seed for generating random intrinsic extinction + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscopic": simulate slitless spectroscopic chips only + # "FGS": simulate FGS chips only (31-42) + # "All": simulate full focal plane + survey_type: "Spectroscopic" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/Volumes/EAGET/C6_data/inputData/" + pointing_file: "pointing_radec_246.5_40.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [80] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [10] + + # Whether to enable astrometric modeling + enable_astrometric_model: False + + # Whether to enable straylight model + enable_straylight_model: True + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + psf_sls_dir: "/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: NO # Whether to add field distortions + add_back: YES # Whether to add sky background + add_dark: YES # Whether to add dark noise + add_readout: YES # Whether to add read-out (Gaussian) noise + add_bias: YES # Whether to add bias-level to images + bias_16channel: YES # Whether to add different biases for 16 channels + gain_16channel: YES # Whether to make different gains for 16 channels + shutter_effect: YES # Whether to add shutter effect + flat_fielding: YES # Whether to add flat-fielding effect + prnu_effect: YES # Whether to add PRNU effect + non_linear: YES # Whether to add non-linearity + cosmic_ray: NO # Whether to add cosmic-ray + cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: YES # Whether to simulate CTE trails + saturbloom: YES # Whether to simulate Saturation & Blooming + add_badcolumns: NO # Whether to add bad columns + add_hotpixels: NO # Whether to add hot pixels + add_deadpixels: NO # Whether to add dead(dark) pixels + bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file