diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index f80123ac920f2892373979f2bc9f39606d1c952c..cb060ab4d65520e8a3ae76b5b5f678a88cfe43d0 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -23,7 +23,7 @@ class Observation(object): self.filter_param = FilterParam() self.Catalog = Catalog - def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None): + def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None, slsPSFOptim=False): # Get WCS for the focal plane if wcs_fp == None: wcs_fp = self.focal_plane.getTanWCS( @@ -34,6 +34,27 @@ class Observation(object): chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) chip.img.wcs = wcs_fp + chip.slsPSFOptim = slsPSFOptim + if chip.chipID in [1, 2, 3, 4, 5, 10, 21, 26, 27, 28, 29, 30] and slsPSFOptim: + chip.img_stack = {} + for id1 in np.arange(2): + gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1] + orders = {} + # for id2 in ['-2','-1','0','1','2']: + for id2 in ['0', '1']: + o_n = "order"+id2 + allbands = {} + for id3 in ['1', '2', '3', '4']: + w_n = "w"+id3 + allbands[w_n] = galsim.ImageF(chip.npix_x, chip.npix_y) + allbands[w_n].setOrigin( + chip.bound.xmin, chip.bound.ymin) + allbands[w_n].wcs = wcs_fp + orders[o_n] = allbands + chip.img_stack[gn] = orders + else: + chip.img_stack = {} + # Get random generators for this chip chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson( seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.) @@ -91,15 +112,17 @@ class Observation(object): input_date_str=date_str, input_time_str=time_str ) - ra_cen, dec_cen = ra_cen[0], dec_cen[0] - ra_offset, dec_offset = pointing.ra - ra_cen, pointing.dec - dec_cen + ra_offset, dec_offset = pointing.ra - \ + ra_cen[0], pointing.dec - dec_cen[0] else: - ra_cen = pointing.ra - dec_cen = pointing.dec ra_offset, dec_offset = 0., 0. + ra_cen = pointing.ra + dec_cen = pointing.dec + slsPSFOpt = False # Prepare necessary chip properties for simulation - chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing) + chip = self.prepare_chip_for_exposure( + chip, ra_cen, dec_cen, pointing, slsPSFOptim=slsPSFOpt) # Initialize SimSteps sim_steps = SimSteps(overall_config=self.config, diff --git a/observation_sim/config/ChipOutput.py b/observation_sim/config/ChipOutput.py index 5efd0616e8870d90a0844279d476d4817a9a2135..7db5e12d465534523c040a7e1d9e66330ece659e 100755 --- a/observation_sim/config/ChipOutput.py +++ b/observation_sim/config/ChipOutput.py @@ -78,12 +78,12 @@ class ChipOutput(object): self.hdr += "\n" self.cat.write(self.hdr) - def cat_add_obj(self, obj, pos_img, pos_shear): + def cat_add_obj(self, obj, pos_img, pos_shear, ra_offset=0., dec_offset=0.): ximg = obj.real_pos.x + 1.0 yimg = obj.real_pos.y + 1.0 line = self.fmt % ( - obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( + obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra + ra_offset, obj.dec + dec_offset, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( self.filt), obj.type, obj.pmra, obj.pmdec, obj.rv, obj.parallax) line += obj.additional_output_str diff --git a/observation_sim/mock_objects/Galaxy.py b/observation_sim/mock_objects/Galaxy.py index c220621e35b1c6121edc2c22b83bb9410e02715e..b0404e7062451b93f4290b7caa2a756bf6d69def 100755 --- a/observation_sim/mock_objects/Galaxy.py +++ b/observation_sim/mock_objects/Galaxy.py @@ -323,28 +323,73 @@ class Galaxy(MockObject): # # if fd_shear is not None: # # gal = gal.shear(fd_shear) - 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)] - starImg.setOrigin(0, 0) + galImg_List = [] + try: + pos_img_local = [0,0] + x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. + y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. + pos_img_local[0] = pos_img.x - x_start + pos_img_local[1] = pos_img.y - y_start + nnx = 0 + nny = 0 + for order in ["A","B"]: + psf, pos_shear = psf_model.get_PSF( + chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos) + star_p = galsim.Convolve(psf, gal) + if nnx == 0: + galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) + nnx = galImg.xmax - galImg.xmin + 1 + nny = galImg.ymax - galImg.ymin + 1 + else: + galImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset) + galImg.setOrigin(0, 0) + # n1 = np.sum(np.isinf(galImg.array)) + # n2 = np.sum(np.isnan(galImg.array)) + # if n1>0 or n2 > 0: + # print("DEBUG: Galaxy, inf:%d, nan:%d"%(n1, n2)) + if np.sum(np.isnan(galImg.array)) > 0: + # ERROR happens + return 2, pos_shear + galImg_List.append(galImg) + for order in ["C","D","E"]: + galImg_List.append(galImg) + except: + psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) + star_p = galsim.Convolve(psf, gal) + galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) + galImg.setOrigin(0, 0) + if np.sum(np.isnan(galImg.array)) > 0: + # ERROR happens + return 2, pos_shear + for order in ["A","B","C","D","E"]: + galImg_List.append(galImg) + + # starImg = gal.drawImage( + # wcs=chip_wcs_local, offset=offset, method='real_space') + + origin_star = [y_nominal - (galImg.center.y - galImg.ymin), + x_nominal - (galImg.center.x - galImg.xmin)] + galImg.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] + gal_end = [origin_star[0] + galImg.array.shape[0] - + 1, origin_star[1] + galImg.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 - subImg_p1 = starImg.array[:, 0:subSlitPos] - star_p1 = galsim.Image(subImg_p1) - star_p1.setOrigin(0, 0) + star_p1s=[] + for galImg in galImg_List: + + subImg_p1 = galImg.array[:, 0:subSlitPos] + star_p1 = galsim.Image(subImg_p1) + star_p1.setOrigin(0, 0) + star_p1s.append(star_p1) origin_p1 = origin_star xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0 ycenter_p1 = y_nominal-0 - sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, + sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -352,21 +397,25 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - # 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) + 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) + + star_p2s=[] + for galImg in galImg_List: - subImg_p2 = starImg.array[:, - subSlitPos+1:starImg.array.shape[1]] - star_p2 = galsim.Image(subImg_p2) - star_p2.setOrigin(0, 0) + subImg_p2 = galImg.array[:, + subSlitPos + 1:galImg.array.shape[1]] + star_p2 = galsim.Image(subImg_p2) + star_p2.setOrigin(0, 0) + star_p2s.append(star_p2) origin_p2 = [origin_star[0], grating_split_pos_chip] xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p2 = y_nominal - 0 - sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, + sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -374,41 +423,41 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - # 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) + 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 elif grating_split_pos_chip <= gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, 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) - 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) + 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, + sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, 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) - 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) + 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 # print(self.y_nominal, starImg.center.y, starImg.ymin) diff --git a/observation_sim/mock_objects/MockObject.py b/observation_sim/mock_objects/MockObject.py index ef660caea1f7a057fff78abe44f1202a9e3745d0..afbbc37a7de9fa2f3e908272a430a1a88d4b497d 100755 --- a/observation_sim/mock_objects/MockObject.py +++ b/observation_sim/mock_objects/MockObject.py @@ -11,6 +11,8 @@ from observation_sim.mock_objects._util import integrate_sed_bandpass, getNormFa getABMAG from observation_sim.mock_objects.SpecDisperser import SpecDisperser +from observation_sim.instruments.chip import chip_utils + class MockObject(object): def __init__(self, param, logger=None): @@ -118,7 +120,6 @@ class MockObject(object): if self.logger: self.logger.error(e) return 2, None - # Set Galsim Parameters if self.getMagFilter(filt) <= 15: folding_threshold = 5.e-4 @@ -195,33 +196,26 @@ class MockObject(object): # DEBUG ######################################################### # print("before convolveGaussXorders, img_s:", img_s) - nan_ids = np.isnan(img_s) - if img_s[nan_ids].shape[0] > 0: - # img_s[nan_ids] = 0 - print("DEBUG: before convolveGaussXorders specImg nan num is", - img_s[nan_ids].shape[0]) + # nan_ids = np.isnan(img_s) + # if img_s[nan_ids].shape[0] > 0: + # # img_s[nan_ids] = 0 + # print("DEBUG: before convolveGaussXorders specImg nan num is", + # img_s[nan_ids].shape[0]) ######################################################### - img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k]) + # img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k]) + orig_off = 0 origin_order_x = v[1] - orig_off origin_order_y = v[2] - orig_off - ######################################################### # DEBUG ######################################################### # 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]) + # 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]) ######################################################### - 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 = galsim.ImageF(img_s) stamp.wcs = local_wcs stamp.setOrigin(origin_order_x, origin_order_y) @@ -230,72 +224,140 @@ class MockObject(object): if bounds.area() == 0: continue chip.img.setOrigin(0, 0) - stamp[bounds] = chip.img[bounds] - chip.sensor.accumulate(photons, stamp) - chip.img[bounds] = stamp[bounds] + chip.img[bounds] = chip.img[bounds]+stamp[bounds] chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) 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 + pos_shear = galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) + if chip.slsPSFOptim: + for k, v in spec_orders.items(): + img_s = v[0] + + 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) + orig_off = [0,0] + origin_order_x = v[1] - orig_off[0] + origin_order_y = v[2] - orig_off[1] + + specImg = galsim.ImageF(img_s) + + 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 + + # orders = {'A': 'order1', 'B': 'order0', 'C': 'order2', 'D': 'order-1', 'E': 'order-2'} + orders = {'A': 'order1', 'B': 'order0', 'C': 'order0', 'D': 'order0', 'E': 'order0'} + gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[1] + if pos_img_local[0] < grating_split_pos: + gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[0] + + + chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(0, 0) + chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] = chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] + specImg[bounds] + chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(chip.bound.xmin, chip.bound.ymin) - # 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 + else: + 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) + + # print('DEBUG: BEGIN -----------',bandNo,k) + + img_s = v[0] + + 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]) + ######################################################### + origin_order_x = v[1] + origin_order_y = v[2] + specImg = galsim.ImageF(img_s) + + specImg.wcs = local_wcs + specImg.setOrigin(origin_order_x, origin_order_y) + try: + specImg = psf_model.get_PSF_AND_convolve_withsubImg(chip, cutImg=specImg, pos_img_local=pos_img_local, bandNo=bandNo, 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 + + 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) + + 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 @@ -359,32 +421,76 @@ class MockObject(object): # 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) - psf_tmp = galsim.Gaussian(sigma=0.002) - star = galsim.Convolve(psf_tmp, star) - starImg = star.drawImage( - nx=60, ny=60, wcs=chip_wcs_local, offset=offset) + # star = galsim.DeltaFunction(gsparams=gsp) + # star = star.withFlux(tel.pupil_area * exptime) + + #psf list :["A","B","C","D","E"] + starImg_List = [] + try: + pos_img_local = [0,0] + x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. + y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. + pos_img_local[0] = pos_img.x - x_start + pos_img_local[1] = pos_img.y - y_start + nnx = 0 + nny = 0 + for order in ["A","B"]: + psf, pos_shear = psf_model.get_PSF( + chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos) + # star_p = galsim.Convolve(psf, star) + star_p = psf.withFlux(tel.pupil_area * exptime) + if nnx == 0: + starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) + nnx = starImg.xmax - starImg.xmin + 1 + nny = starImg.ymax - starImg.ymin + 1 + else: + starImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset) + # n1 = np.sum(np.isinf(starImg.array)) + # n2 = np.sum(np.isnan(starImg.array)) + # if n1>0 or n2 > 0: + # print("DEBUG: MockObject, inf:%d, nan:%d"%(n1, n2)) + starImg.setOrigin(0, 0) + starImg_List.append(starImg) + for order in ["C","D","E"]: + starImg_List.append(starImg) + except: + psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) + # star_p = galsim.Convolve(psf, star) + star_p = psf.withFlux(tel.pupil_area * exptime) + starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) + starImg.setOrigin(0, 0) + for order in ["A","B","C","D","E"]: + starImg_List.append(starImg) + + + # psf_tmp = galsim.Gaussian(sigma=0.002) + # star = galsim.Convolve(psf_tmp, star) + + # 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 + star_p1s=[] + for starImg in starImg_List: - subImg_p1 = starImg.array[:, 0:subSlitPos] - star_p1 = galsim.Image(subImg_p1) + subImg_p1 = starImg.array[:, 0:subSlitPos] + star_p1 = galsim.Image(subImg_p1) + star_p1.setOrigin(0, 0) + star_p1s.append(star_p1) origin_p1 = origin_star - star_p1.setOrigin(0, 0) xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p1 = y_nominal - 0 - sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, + sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -392,20 +498,25 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - # 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) + 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) + + + star_p2s=[] + for starImg in starImg_List: - subImg_p2 = starImg.array[:, + subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]] - star_p2 = galsim.Image(subImg_p2) - star_p2.setOrigin(0, 0) + star_p2 = galsim.Image(subImg_p2) + star_p2.setOrigin(0, 0) + star_p2s.append(star_p2) origin_p2 = [origin_star[0], grating_split_pos_chip] xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p2 = y_nominal - 0 - sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, + sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -413,38 +524,38 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - # 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) + 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 elif grating_split_pos_chip <= gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, 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) - 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) + 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, + sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, 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) - 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) + 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 return 1, pos_shear diff --git a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py index 6a1e3ef05c70fa418db91149eb1e60bf9dc8b81a..2b2d7c29baa9e9965725da0cf92fac98b8495ae2 100644 --- a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py +++ b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py @@ -65,10 +65,30 @@ class SpecDisperser(object): # self.img_x = orig_img.shape[1] # self.img_y = orig_img.shape[0] - self.thumb_img = np.abs(orig_img.array) - self.thumb_x = orig_img.center.x - self.thumb_y = orig_img.center.y - self.img_sh = orig_img.array.shape + # 5 orders, A, B , + orderName=["A","B","C","D","E"] + self.orig_img_orders = OrderedDict() + if isinstance(orig_img, list): + orig_img_list = orig_img + list_len = len(orig_img_list) + if list_len < 5: + for i in np.arange(5-list_len): + orig_img_list.append(orig_img_list[list_len-1]) + for i, k in enumerate(orig_img_list): + self.orig_img_orders[orderName[i]] = k + + + if isinstance(orig_img, galsim.Image): + for i in np.arange(5): + self.orig_img_orders[orderName[i]] = orig_img + + + orig_img_one = self.orig_img_orders["A"] + + self.thumb_img = np.abs(orig_img_one.array) + self.thumb_x = orig_img_one.center.x + self.thumb_y = orig_img_one.center.y + self.img_sh = orig_img_one.array.shape self.id = gid @@ -78,10 +98,13 @@ class SpecDisperser(object): self.isAlongY = isAlongY self.flat_cube = flat_cube if self.isAlongY == 1: - self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img.center.x, - yc=orig_img.center.y, isClockwise=1) + for order in orderName: + self.orig_img_orders[order], self.thumb_x, self.thumb_y = rotate90(array_orig=self.orig_img_orders[order], xc=orig_img_one.center.x, + yc=orig_img_one.center.y, isClockwise=1) + # self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img_one.center.x, + # yc=orig_img_one.center.y, isClockwise=1) - self.img_sh = orig_img.array.T.shape + self.img_sh = self.orig_img_orders[order].array.T.shape self.xcenter = ycenter self.ycenter = xcenter @@ -111,10 +134,16 @@ class SpecDisperser(object): def compute_spec(self, beam): + # if beam == "B": + # return self.thumb_img, self.origin[1], self.origin[0], None, None, None + from .disperse_c import interp from .disperse_c import disperse # from MockObject.disperse_c import disperse - + self.thumb_img = np.abs(self.orig_img_orders[beam].array) + self.thumb_x = self.orig_img_orders[beam].center.x + self.thumb_y = self.orig_img_orders[beam].center.y + self.img_sh = self.orig_img_orders[beam].array.shape dx = self.grating_conf.dxlam[beam] xoff = 0 ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff), @@ -169,7 +198,8 @@ class SpecDisperser(object): dyc = cast[int](np.floor(ytrace_beam+0.5)) - dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5)) + # dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5)) + dypix = dyc - dyc[0] + x0[0] frac_ids = yfrac_beam < 0 @@ -248,7 +278,8 @@ class SpecDisperser(object): # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), - flat_index[nonz], yfrac_beam[nonz], + flat_index[nonz], + yfrac_beam[nonz], sensitivity_beam[nonz], modelf, x0, array(self.img_sh, @@ -258,11 +289,24 @@ class SpecDisperser(object): lam_beam[lam_index][nonz]) model = modelf.reshape(beam_sh) + # n1 = np.sum(np.isinf(model)) + # n2 = np.sum(np.isnan(model)) + # n3 = np.sum(np.isinf(modelf)) + # n4 = np.sum(np.isnan(modelf)) + # if n1>0 or n2 > 0: + # print("DEBUG: SpecDisperser, inf:%d, nan:%d--------%d,%d"%(n1, n2, n3, n4)) + # print(dypix) + # n1 = np.sum(np.isinf(self.thumb_img.astype(np.float32))) + # n2 = np.sum(np.isnan(self.thumb_img.astype(np.float32))) + # n3 = np.sum(np.isinf(yfrac_beam)) + # n4 = np.sum(np.isnan(yfrac_beam)) + # n5 = np.sum(np.isinf(sensitivity_beam)) + # n6 = np.sum(np.isnan(sensitivity_beam)) + # print("DEBUG: SpecDisperser, innput ---inf:%d, nan:%d, yfrac_beam:%d/%d, sensitivity_beam:%d/%d"%(n1, n2, n3, n4, n5, n6)) self.beam_flux[beam] = sum(modelf) if self.isAlongY == 1: model, _, _ = rotate90(array_orig=model, isClockwise=0) - return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None): diff --git a/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx b/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx index 41dbf1402e7aa4d756865489898cc738187a094a..707c6707fa2f3044f915bf54c6f8087970452ee2 100644 --- a/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx +++ b/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx @@ -20,7 +20,30 @@ import cython cdef extern from "math.h": double sqrt(double x) double exp(double x) + + +def check_nan2D(np.ndarray[FTYPE_t, ndim=2] arr): + cdef int i, j + cdef int nrows = arr.shape[0] + cdef int ncols = arr.shape[1] + # 遍历数组的每个元素并检查是否存在 NaN + for i in range(nrows): + for j in range(ncols): + if np.isnan(arr[i, j]) | np.isinf(arr[i, j]): + return True + return False + +def check_nan1d(np.ndarray[DTYPE_t, ndim=1] arr): + cdef int i + cdef int n = arr.shape[0] + + # 遍历数组的每个元素并检查是否存在 NaN + for i in range(n): + if np.isnan(arr[i]) | np.isinf(arr[i]): + return True + return False + @cython.boundscheck(False) @cython.wraparound(False) @cython.embedsignature(True) @@ -53,6 +76,18 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, nk = len(idxl) nl = len(full) + + + #if check_nan2D(flam): + # print("DEBUG: disperse, input Array 'flam' contains NaN.") + #if check_nan1d(ysens): + # print("DEBUG: disperse, input Array 'ysens' contains NaN.") + #if check_nan1d(yfrac): + # print("DEBUG: disperse, input Array 'yfrac' contains NaN.") + #if check_nan1d(full): + # print("DEBUG: disperse, input Array 'full' contains NaN before processing.") + + if (flat is not None): nlamb = len(wlambda) @@ -95,14 +130,15 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, else: for i in range(0-x0[1], x0[1]): - if (x0[1]+i < 0) | (x0[1]+i >= shd[1]): + x_pos = x0[1]+i + if (x_pos < 0) | (x_pos >= shd[1]): continue for j in range(0-x0[0], x0[0]): - if (x0[0]+j < 0) | (x0[0]+j >= shd[0]): + y_pos = x0[0]+j + if (y_pos < 0) | (y_pos >= shd[0]): continue - - fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17 + fl_ij = flam[y_pos, x_pos] #/1.e-17 if (fl_ij == 0): continue @@ -110,10 +146,13 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, k1 = idxl[k]+j*shg[1]+i if (k1 >= 0) & (k1 < nl): full[k1] += ysens[k]*fl_ij*(1-yfrac[k]) - k2 = idxl[k]+(j+1)*shg[1]+i if (k2 >= 0) & (k2 < nl): full[k2] += ysens[k]*fl_ij*yfrac[k] + + #if (check_nan1d(full)): + # print("DEBUG: disperse, output Array 'full' contains NaN after processing.+++++++++++++++++++++++++++") + return True diff --git a/observation_sim/psf/PSFGauss.py b/observation_sim/psf/PSFGauss.py index a21b0c8963da1aa7b4e22701d06b234b8402a69c..e3036d9293c37397505a4abf65b367b5d9429a8f 100644 --- a/observation_sim/psf/PSFGauss.py +++ b/observation_sim/psf/PSFGauss.py @@ -17,7 +17,7 @@ class PSFGauss(PSFModel): self.fwhm = self.fwhmGauss(r=psfRa) self.sigGauss = psfRa # 80% light radius self.sigSpin = sigSpin - self.psf = galsim.Gaussian(flux=1.0, fwhm=fwhm) + self.psf = galsim.Gaussian(flux=1.0, fwhm=self.fwhm) def perfGauss(self, r, sig): """ @@ -122,4 +122,4 @@ class PSFGauss(PSFModel): # return ell, beta, qr PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians) - return self.psf.shear(PSFshear), PSFshear + return self.psf.shear(PSFshear), PSFshear \ No newline at end of file diff --git a/observation_sim/psf/PSFInterpSLS.py b/observation_sim/psf/PSFInterpSLS.py index 7adfedc5b953339492940d41b08e0f31c062213c..90ac027a4ce4d38a5fc93298d717b3b197e89351 100644 --- a/observation_sim/psf/PSFInterpSLS.py +++ b/observation_sim/psf/PSFInterpSLS.py @@ -20,8 +20,10 @@ import os from astropy.io import fits from astropy.modeling.models import Gaussian2D -from scipy import signal - +from scipy import signal, interpolate +import datetime +import gc +# from jax import numpy as jnp LOG_DEBUG = False # ***# NPSF = 900 # ***# 30*30 @@ -433,7 +435,16 @@ class PSFInterpSLS(PSFModel): # 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) + ###DEBGU + ids_szero = PSF_int_trans<0 + n01 = PSF_int_trans[ids_szero].shape[0] + + n1 = np.sum(np.isinf(PSF_int_trans)) + n2 = np.sum(np.isnan(PSF_int_trans)) + if n1>0 or n2>0: + print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d"%(n1, n2, n01)) + #### # from astropy.io import fits # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans) @@ -479,6 +490,215 @@ class PSFInterpSLS(PSFModel): return PSF_int_trans, PSF_int + def get_PSF_AND_convolve_withsubImg(self, chip, cutImg=None, pos_img_local=[1000, 1000], bandNo=1, 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) + # centerPos_local = cutImg.ncol/2. + 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 + pos_p = psf_b_dat[1].data/chip.pix_size - np.array([y_start, x_start]) + + pc_coeff = psf_b_dat[2].data + pcs = psf_b_dat[0].data + + npc = 10 + m_size = int(pcs.shape[0]**0.5) + sumImg = np.sum(cutImg.array) + tmp_img = cutImg*0 + for j in np.arange(npc): + X_ = np.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) + Z_ = (pc_coeff[j].astype(np.float32)).flatten() + # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) + cx_len = int(chip.npix_x) + cy_len = int(chip.npix_y) + n_x = np.arange(0, cx_len, 1, dtype = int) + n_y = np.arange(0, cy_len, 1, dtype = int) + M, N = np.meshgrid(n_x, n_y) + # t1=datetime.datetime.now() + # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), + # method='nearest',fill_value=1.0) + b_img = galsim.Image(cx_len, cy_len) + b_img.setOrigin(0,0) + bounds = cutImg.bounds & b_img.bounds + if bounds.area() == 0: + continue + + # ys = cutImg.ymin + # if ys < 0: + # ys = 0 + # ye = cutImg.ymin+cutImg.nrow + # if ye >= cy_len-1: + # ye = cy_len-1 + # if ye - ys <=0: + # continue + # xs = cutImg.xmin + # if xs < 0: + # xs = 0 + # xe = cutImg.xmin+cutImg.ncol + # if xe >= cx_len-1: + # xe = cx_len-1 + # if xe - xs <=0: + # continue + ys = bounds.ymin + ye = bounds.ymax+1 + xs = bounds.xmin + xe = bounds.xmax+1 + U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe],N[ys:ye, xs:xe]), + method='nearest',fill_value=1.0) + # t2=datetime.datetime.now() + + # print("time interpolate:", t2-t1) + + # if U.shape != cutImg.array.shape: + # print('DEBUG:SHAPE',cutImg.ncol,cutImg.nrow,cutImg.xmin, cutImg.ymin) + # continue + img_tmp = cutImg + img_tmp[bounds] = img_tmp[bounds]*U + psf = pcs[:, j].reshape(m_size, m_size) + tmp_img = tmp_img + signal.fftconvolve(img_tmp.array, psf, mode='same', axes=None) + + # t3=datetime.datetime.now() + # print("time convole:", t3-t2) + del U + del img_tmp + if np.sum(tmp_img.array)==0: + tmp_img = cutImg + else: + tmp_img = tmp_img/np.sum(tmp_img.array)*sumImg + return tmp_img + + + def convolveFullImgWithPCAPSF(self, chip, folding_threshold=5.e-3): + keys_L1= chip_utils.getChipSLSGratingID(chip.chipID) + # keys_L2 = ['order-2','order-1','order0','order1','order2'] + keys_L2 = ['order0','order1'] + keys_L3 = ['w1','w2','w3','w4'] + + npca = 10 + + x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. + y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. + + for i,gt in enumerate(keys_L1): + psfCo = self.grating1_data + if i > 0: + psfCo = self.grating2_data + for od in keys_L2: + psfCo_L2 = psfCo['order1'] + if od in ['order-2','order-1','order0','order2']: + psfCo_L2 = psfCo['order0'] + for w in keys_L3: + img = chip.img_stack[gt][od][w] + pcs = psfCo_L2['band'+w[1]]['band_data'][0].data + pos_p = psfCo_L2['band'+w[1]]['band_data'][1].data/chip.pix_size - np.array([y_start, x_start]) + pc_coeff = psfCo_L2['band'+w[1]]['band_data'][2].data + # print("DEBUG-----------",np.max(pos_p[:,1]),np.min(pos_p[:,1]), np.max(pos_p[:,0]),np.min(pos_p[:,0])) + sum_img = np.sum(img.array) + + + # coeff_mat = np.zeros([npca, chip.npix_y, chip.npix_x]) + # for m in np.arange(chip.npix_y): + # for n in np.arange(chip.npix_x): + # px = n + # py = m + + # 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, 3]) + # pc_coeff_4p = np.zeros([npca, 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] + # # print(pos_p[smaller_ids, 1],pos_p[smaller_ids, 0]) + # nearest4p[i, 2] = dits2_sortlist[i][1] + # pc_coeff_4p[:, i] = pc_coeff[npca, smaller_ids] + # # idw_dist = 1/(np.sqrt((px-nearest4p[:, 0]) * (px-nearest4p[:, 0]) + ( + # # py-nearest4p[:, 1]) * (py-nearest4p[:, 1]))) + # idw_dist = 1/(np.sqrt(nearest4p[:, 2])) + + # coeff_int = np.zeros(npca) + # for i in np.arange(4): + # coeff_int = coeff_int + pc_coeff_4p[:, i]*idw_dist[i] + # coeff_mat[:, m, n] = coeff_int + + m_size = int(pcs.shape[0]**0.5) + tmp_img = np.zeros_like(img.array,dtype=np.float32) + for j in np.arange(npca): + print(gt, od, w, j) + X_ = np.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) + Z_ = (pc_coeff[j].astype(np.float32)).flatten() + # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) + sub_size = 4 + cx_len = int(chip.npix_x/sub_size) + cy_len = int(chip.npix_y/sub_size) + n_x = np.arange(0, chip.npix_x, sub_size, dtype = int) + n_y = np.arange(0, chip.npix_y, sub_size, dtype = int) + + M, N = np.meshgrid(n_x, n_y) + t1=datetime.datetime.now() + # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), + # method='nearest',fill_value=1.0) + U1 = interpolate.griddata(X_, Z_, (M, N), + method='nearest',fill_value=1.0) + U = np.zeros_like(chip.img.array, dtype=np.float32) + for mi in np.arange(cy_len): + for mj in np.arange(cx_len): + U[mi*sub_size:(mi+1)*sub_size, mj*sub_size:(mj+1)*sub_size]=U1[mi,mj] + t2=datetime.datetime.now() + + print("time interpolate:", t2-t1) + + img_tmp = img.array*U + psf = pcs[:, j].reshape(m_size, m_size) + tmp_img = tmp_img + signal.fftconvolve(img_tmp, psf, mode='same', axes=None) + + t3=datetime.datetime.now() + print("time convole:", t3-t2) + del U + del U1 + + chip.img = chip.img + tmp_img*sum_img/np.sum(tmp_img) + del tmp_img + gc.collect() + # 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' diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 741d9639c6cabb31626884eae3959348fa40748e..6c9c7c197577b2c8a3fda456043366b03d40691e 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -194,7 +194,8 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): if isUpdated == 1: # TODO: add up stats - self.chip_output.cat_add_obj(obj, pos_img, pos_shear) + self.chip_output.cat_add_obj( + obj, pos_img, pos_shear, ra_offset=self.ra_offset, dec_offset=self.dec_offset) pass elif isUpdated == 0: missed_obj += 1 @@ -217,6 +218,26 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): obj.unload_SED() del obj # gc.collect() + + if chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"] and chip.slsPSFOptim: + # from observation_sim.instruments.chip import chip_utils as chip_utils + # gn = chip_utils.getChipSLSGratingID(chip.chipID)[0] + # img1 = np.zeros([2,chip.img.array.shape[0],chip.img.array.shape[1]]) + + # for id1 in np.arange(2): + # gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1] + # img_i = 0 + # for id2 in ['0','1']: + # o_n = "order"+id2 + # for id3 in ['1','2','3','4']: + # w_n = "w"+id3 + # img1[img_i] = img1[img_i] + chip.img_stack[gn][o_n][w_n].array + # img_i = img_i + 1 + # from astropy.io import fits + # fits.writeto('order0.fits',img1[0],overwrite=True) + # fits.writeto('order1.fits',img1[1],overwrite=True) + + psf_model.convolveFullImgWithPCAPSF(chip) del psf_model gc.collect()