Skip to content
GitLab
Commit 1ca08d23 authored by Zhang Xin's avatar Zhang Xin
Browse files

add sls psf model

parent 278305de
Hide whitespace changes
Inline Side-by-side
ObservationSim/Config/Config.py
View file @ 1ca08d23
...@@ -20,6 +20,7 @@ def config_dir(config, work_dir=None, data_dir=None): ...@@ -20,6 +20,7 @@ def config_dir(config, work_dir=None, data_dir=None):
# PSF data directory # PSF data directory
if config["psf_setting"]["psf_model"] == "Interp": 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_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 return path_dict
......
ObservationSim/Instrument/Chip/ChipUtils.py
View file @ 1ca08d23
...@@ -21,6 +21,22 @@ def log_info(msg, logger=None): ...@@ -21,6 +21,22 @@ def log_info(msg, logger=None):
else: else:
print(msg, flush=True) 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): def getChipSLSConf(chipID):
confFile = '' confFile = ''
if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf'] if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
......
ObservationSim/MockObject/Galaxy.py
View file @ 1ca08d23
...@@ -248,10 +248,27 @@ class Galaxy(MockObject): ...@@ -248,10 +248,27 @@ class Galaxy(MockObject):
xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos 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)): 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 = 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_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape) disk = disk.shear(disk_shape)
...@@ -272,14 +289,14 @@ class Galaxy(MockObject): ...@@ -272,14 +289,14 @@ class Galaxy(MockObject):
g2 += fd_shear.g2 g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2) gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear) 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 # if not big_galaxy: # Not apply PSF for very big galaxy
gal = galsim.Convolve(psf, gal) # gal = galsim.Convolve(psf, gal)
# if fd_shear is not None: # # if fd_shear is not None:
# gal = gal.shear(fd_shear) # # 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), origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)] x_nominal - (starImg.center.x - starImg.xmin)]
...@@ -301,12 +318,16 @@ class Galaxy(MockObject): ...@@ -301,12 +318,16 @@ class Galaxy(MockObject):
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, 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], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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]] subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
...@@ -318,12 +339,16 @@ class Galaxy(MockObject): ...@@ -318,12 +339,16 @@ class Galaxy(MockObject):
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, 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], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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_p1
del sdp_p2 del sdp_p2
...@@ -331,25 +356,33 @@ class Galaxy(MockObject): ...@@ -331,25 +356,33 @@ class Galaxy(MockObject):
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, 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], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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 del sdp
elif grating_split_pos_chip>=gal_end[1]: elif grating_split_pos_chip>=gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, 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], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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 del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin) # print(self.y_nominal, starImg.center.y, starImg.ymin)
del psf # del psf
return 1, pos_shear return 1, pos_shear
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.): def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
......
ObservationSim/MockObject/MockObject.py
View file @ 1ca08d23
...@@ -4,7 +4,7 @@ import astropy.constants as cons ...@@ -4,7 +4,7 @@ import astropy.constants as cons
from astropy import wcs from astropy import wcs
from astropy.table import Table 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, \ from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, \
getABMAG getABMAG
from ObservationSim.MockObject.SpecDisperser import SpecDisperser from ObservationSim.MockObject.SpecDisperser import SpecDisperser
...@@ -222,9 +222,69 @@ class MockObject(object): ...@@ -222,9 +222,69 @@ class MockObject(object):
del stamp del stamp
del spec_orders 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, 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): exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
if normFilter is not None: if normFilter is not None:
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed, sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
...@@ -262,22 +322,38 @@ class MockObject(object): ...@@ -262,22 +322,38 @@ class MockObject(object):
xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062, 'C': 4.035447379743442, xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062, 'C': 4.035447379743442,
'D': 5.5684364343742825, 'E': 16.260021029735388} 'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos 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)): for i in range(len(bandpass_list)):
bandpass = bandpass_list[i] # bandpass = bandpass_list[i]
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, brange = branges[i]
folding_threshold=folding_threshold)
# 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 = galsim.DeltaFunction(gsparams=gsp)
star = star.withFlux(tel.pupil_area * exptime) 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), origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)] x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0,0) starImg.setOrigin(0,0)
gal_origin = [origin_star[0], origin_star[1]] 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] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]: if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse ## part img disperse
...@@ -292,12 +368,15 @@ class MockObject(object): ...@@ -292,12 +368,15 @@ class MockObject(object):
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, 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], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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]] subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
...@@ -309,12 +388,15 @@ class MockObject(object): ...@@ -309,12 +388,15 @@ class MockObject(object):
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, 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], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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_p1
del sdp_p2 del sdp_p2
...@@ -322,23 +404,29 @@ class MockObject(object): ...@@ -322,23 +404,29 @@ class MockObject(object):
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, 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], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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 sdp
elif grating_split_pos_chip >= gal_end[1]: elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, 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], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) 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 sdp
del psf # del psf
return 1, pos_shear return 1, pos_shear
def SNRestimate(self, img_obj, flux, noise_level=0.0, seed=31415): def SNRestimate(self, img_obj, flux, noise_level=0.0, seed=31415):
......
ObservationSim/MockObject/SpecDisperser/SpecDisperser.py
View file @ 1ca08d23
...@@ -155,7 +155,7 @@ class SpecDisperser(object): ...@@ -155,7 +155,7 @@ class SpecDisperser(object):
sensitivity_beam = ysens sensitivity_beam = ysens
len_spec_x = len(dx) 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) beam_sh = (self.img_sh[0] + len_spec_y, self.img_sh[1] + len_spec_x)
modelf = zeros(product(beam_sh), dtype=float) modelf = zeros(product(beam_sh), dtype=float)
......
ObservationSim/MockObject/_util.py
View file @ 1ca08d23
...@@ -571,4 +571,14 @@ def convolveGaussXorders(img=None, sigma = 1): ...@@ -571,4 +571,14 @@ def convolveGaussXorders(img=None, sigma = 1):
convImg = signal.fftconvolve(img, psf, mode='full', axes=None) convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
return convImg, offset 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
ObservationSim/ObservationSim.py
View file @ 1ca08d23
...@@ -15,7 +15,7 @@ from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensio ...@@ -15,7 +15,7 @@ from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensio
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from ObservationSim.Instrument.Chip import Effects from ObservationSim.Instrument.Chip import Effects
from ObservationSim.Straylight import calculateSkyMap_split_g 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._util import get_shear_field, makeSubDir_PointingList
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
...@@ -63,7 +63,10 @@ class Observation(object): ...@@ -63,7 +63,10 @@ class Observation(object):
if self.config["psf_setting"]["psf_model"] == "Gauss": if self.config["psf_setting"]["psf_model"] == "Gauss":
psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"])
elif self.config["psf_setting"]["psf_model"] == "Interp": 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: else:
chip_output.Log_error("unrecognized PSF model type!!", flush=True) chip_output.Log_error("unrecognized PSF model type!!", flush=True)
...@@ -198,6 +201,10 @@ class Observation(object): ...@@ -198,6 +201,10 @@ class Observation(object):
obj = self.cat.objs[j] 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 # load and convert SED; also caculate object's magnitude in all CSST bands
try: try:
sed_data = self.cat.load_sed(obj) sed_data = self.cat.load_sed(obj)
......
ObservationSim/PSF/PSFInterpSLS.py 0 → 100644
View file @ 1ca08d23
'''
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')
ObservationSim/PSF/__init__.py
View file @ 1ca08d23
...@@ -2,4 +2,5 @@ from .PSFModel import PSFModel ...@@ -2,4 +2,5 @@ from .PSFModel import PSFModel
from .PSFGauss import PSFGauss from .PSFGauss import PSFGauss
# from .PSFInterp.PSFInterp import PSFInterp # from .PSFInterp.PSFInterp import PSFInterp
from .PSFInterp import PSFInterp from .PSFInterp import PSFInterp
from .PSFInterpSLS import PSFInterpSLS
from .FieldDistortion import FieldDistortion from .FieldDistortion import FieldDistortion
\ No newline at end of file
config/config_C6_dev.yaml 0 → 100644
View file @ 1ca08d23
---
###############################################
#
# 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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment