Skip to content
GitLab
Commit 398f6272 authored by Zhang Xin's avatar Zhang Xin
Browse files

modify sls psf convolve mode, first convolve psf, then pixel to image; fix Gaus psf bug

parent 45494b69
Hide whitespace changes
Inline Side-by-side
observation_sim/ObservationSim.py
View file @ 398f6272
...@@ -116,7 +116,7 @@ class Observation(object): ...@@ -116,7 +116,7 @@ class Observation(object):
ra_cen = pointing.ra ra_cen = pointing.ra
dec_cen = pointing.dec dec_cen = pointing.dec
slsPSFOpt = True slsPSFOpt = False
# Prepare necessary chip properties for simulation # Prepare necessary chip properties for simulation
chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing, slsPSFOptim = slsPSFOpt) chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing, slsPSFOptim = slsPSFOpt)
......
observation_sim/mock_objects/Galaxy.py
View file @ 398f6272
...@@ -323,28 +323,73 @@ class Galaxy(MockObject): ...@@ -323,28 +323,73 @@ class Galaxy(MockObject):
# # 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( galImg_List = []
wcs=chip_wcs_local, offset=offset, method='real_space') try:
pos_img_local = [0,0]
origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
x_nominal - (starImg.center.x - starImg.xmin)] y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
starImg.setOrigin(0, 0) 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_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - gal_end = [origin_star[0] + galImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1] 1, origin_star[1] + galImg.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
subImg_p1 = starImg.array[:, 0:subSlitPos] star_p1s=[]
star_p1 = galsim.Image(subImg_p1) for galImg in galImg_List:
star_p1.setOrigin(0, 0)
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 origin_p1 = origin_star
xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0 xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-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, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -352,21 +397,25 @@ class Galaxy(MockObject): ...@@ -352,21 +397,25 @@ class Galaxy(MockObject):
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], # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
psf_model=psf_model, bandNo=i + 1, # psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, # grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # local_wcs=chip_wcs_local, pos_img=pos_img)
star_p2s=[]
for galImg in galImg_List:
subImg_p2 = starImg.array[:, subImg_p2 = galImg.array[:,
subSlitPos+1:starImg.array.shape[1]] subSlitPos + 1:galImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0) star_p2.setOrigin(0, 0)
star_p2s.append(star_p2)
origin_p2 = [origin_star[0], grating_split_pos_chip] origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 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, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -374,41 +423,41 @@ class Galaxy(MockObject): ...@@ -374,41 +423,41 @@ class Galaxy(MockObject):
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], # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1, # psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, # grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1 del sdp_p1
del sdp_p2 del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]: 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, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], 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], # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, # psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, # grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # 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=galImg_List, 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=brange[0], band_end=brange[1], 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], # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, # psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, # grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin) # print(self.y_nominal, starImg.center.y, starImg.ymin)
......
observation_sim/mock_objects/MockObject.py
View file @ 398f6272
...@@ -196,33 +196,26 @@ class MockObject(object): ...@@ -196,33 +196,26 @@ class MockObject(object):
# DEBUG # DEBUG
######################################################### #########################################################
# print("before convolveGaussXorders, img_s:", img_s) # print("before convolveGaussXorders, img_s:", img_s)
nan_ids = np.isnan(img_s) # nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0: # if img_s[nan_ids].shape[0] > 0:
# img_s[nan_ids] = 0 # # img_s[nan_ids] = 0
print("DEBUG: before convolveGaussXorders specImg nan num is", # print("DEBUG: before convolveGaussXorders specImg nan num is",
img_s[nan_ids].shape[0]) # 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_x = v[1] - orig_off
origin_order_y = v[2] - orig_off origin_order_y = v[2] - orig_off
######################################################### #########################################################
# DEBUG # DEBUG
######################################################### #########################################################
# print("DEBUG: orig_off is", orig_off) # print("DEBUG: orig_off is", orig_off)
nan_ids = np.isnan(img_s) # nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0: # if img_s[nan_ids].shape[0] > 0:
img_s[nan_ids] = 0 # img_s[nan_ids] = 0
print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0]) # print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
######################################################### #########################################################
specImg = galsim.ImageF(img_s) stamp = 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.wcs = local_wcs
stamp.setOrigin(origin_order_x, origin_order_y) stamp.setOrigin(origin_order_x, origin_order_y)
...@@ -231,9 +224,7 @@ class MockObject(object): ...@@ -231,9 +224,7 @@ class MockObject(object):
if bounds.area() == 0: if bounds.area() == 0:
continue continue
chip.img.setOrigin(0, 0) chip.img.setOrigin(0, 0)
stamp[bounds] = chip.img[bounds] chip.img[bounds] = chip.img[bounds]+stamp[bounds]
chip.sensor.accumulate(photons, stamp)
chip.img[bounds] = stamp[bounds]
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp del stamp
del spec_orders del spec_orders
...@@ -430,32 +421,76 @@ class MockObject(object): ...@@ -430,32 +421,76 @@ class MockObject(object):
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
# folding_threshold=folding_threshold) # 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( # star = galsim.DeltaFunction(gsparams=gsp)
nx=60, ny=60, wcs=chip_wcs_local, offset=offset) # 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), 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)
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] - gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1] 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
star_p1s=[]
for starImg in starImg_List:
subImg_p1 = starImg.array[:, 0:subSlitPos] subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1) star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0)
star_p1s.append(star_p1)
origin_p1 = origin_star origin_p1 = origin_star
star_p1.setOrigin(0, 0)
xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p1 = y_nominal - 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, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -463,20 +498,25 @@ class MockObject(object): ...@@ -463,20 +498,25 @@ class MockObject(object):
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], # 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, # psf_model=psf_model, bandNo=i+1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # 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]] subSlitPos + 1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0) star_p2.setOrigin(0, 0)
star_p2s.append(star_p2)
origin_p2 = [origin_star[0], grating_split_pos_chip] origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 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, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -484,38 +524,38 @@ class MockObject(object): ...@@ -484,38 +524,38 @@ class MockObject(object):
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], # 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, # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1 del sdp_p1
del sdp_p2 del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]: 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, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], 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], # 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, # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # 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_List, 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=brange[0], band_end=brange[1], 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], # 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, # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img) # 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
......
observation_sim/mock_objects/SpecDisperser/SpecDisperser.py
View file @ 398f6272
...@@ -65,10 +65,30 @@ class SpecDisperser(object): ...@@ -65,10 +65,30 @@ class SpecDisperser(object):
# self.img_x = orig_img.shape[1] # self.img_x = orig_img.shape[1]
# self.img_y = orig_img.shape[0] # self.img_y = orig_img.shape[0]
self.thumb_img = np.abs(orig_img.array) # 5 orders, A, B ,
self.thumb_x = orig_img.center.x orderName=["A","B","C","D","E"]
self.thumb_y = orig_img.center.y self.orig_img_orders = OrderedDict()
self.img_sh = orig_img.array.shape 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 self.id = gid
...@@ -78,10 +98,13 @@ class SpecDisperser(object): ...@@ -78,10 +98,13 @@ class SpecDisperser(object):
self.isAlongY = isAlongY self.isAlongY = isAlongY
self.flat_cube = flat_cube self.flat_cube = flat_cube
if self.isAlongY == 1: if self.isAlongY == 1:
self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img.center.x, for order in orderName:
yc=orig_img.center.y, isClockwise=1) 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.xcenter = ycenter
self.ycenter = xcenter self.ycenter = xcenter
...@@ -111,10 +134,16 @@ class SpecDisperser(object): ...@@ -111,10 +134,16 @@ class SpecDisperser(object):
def compute_spec(self, beam): 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 interp
from .disperse_c import disperse from .disperse_c import disperse
# from MockObject.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] dx = self.grating_conf.dxlam[beam]
xoff = 0 xoff = 0
ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff), 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): ...@@ -169,7 +198,8 @@ class SpecDisperser(object):
dyc = cast[int](np.floor(ytrace_beam+0.5)) 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 frac_ids = yfrac_beam < 0
...@@ -248,7 +278,8 @@ class SpecDisperser(object): ...@@ -248,7 +278,8 @@ class SpecDisperser(object):
# beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]
status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), 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], sensitivity_beam[nonz],
modelf, x0, modelf, x0,
array(self.img_sh, array(self.img_sh,
...@@ -258,11 +289,24 @@ class SpecDisperser(object): ...@@ -258,11 +289,24 @@ class SpecDisperser(object):
lam_beam[lam_index][nonz]) lam_beam[lam_index][nonz])
model = modelf.reshape(beam_sh) 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) self.beam_flux[beam] = sum(modelf)
if self.isAlongY == 1: if self.isAlongY == 1:
model, _, _ = rotate90(array_orig=model, isClockwise=0) model, _, _ = rotate90(array_orig=model, isClockwise=0)
return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens
def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None): def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None):
......
observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx
View file @ 398f6272
...@@ -20,7 +20,30 @@ import cython ...@@ -20,7 +20,30 @@ import cython
cdef extern from "math.h": cdef extern from "math.h":
double sqrt(double x) double sqrt(double x)
double exp(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.boundscheck(False)
@cython.wraparound(False) @cython.wraparound(False)
@cython.embedsignature(True) @cython.embedsignature(True)
...@@ -53,6 +76,18 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -53,6 +76,18 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
nk = len(idxl) nk = len(idxl)
nl = len(full) 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): if (flat is not None):
nlamb = len(wlambda) nlamb = len(wlambda)
...@@ -95,14 +130,15 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -95,14 +130,15 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
else: else:
for i in range(0-x0[1], x0[1]): 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 continue
for j in range(0-x0[0], x0[0]): 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 continue
fl_ij = flam[y_pos, x_pos] #/1.e-17
fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17
if (fl_ij == 0): if (fl_ij == 0):
continue continue
...@@ -110,10 +146,13 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -110,10 +146,13 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
k1 = idxl[k]+j*shg[1]+i k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl): if (k1 >= 0) & (k1 < nl):
full[k1] += ysens[k]*fl_ij*(1-yfrac[k]) full[k1] += ysens[k]*fl_ij*(1-yfrac[k])
k2 = idxl[k]+(j+1)*shg[1]+i k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl): if (k2 >= 0) & (k2 < nl):
full[k2] += ysens[k]*fl_ij*yfrac[k] full[k2] += ysens[k]*fl_ij*yfrac[k]
#if (check_nan1d(full)):
# print("DEBUG: disperse, output Array 'full' contains NaN after processing.+++++++++++++++++++++++++++")
return True return True
......
observation_sim/psf/PSFGauss.py
View file @ 398f6272
...@@ -17,7 +17,7 @@ class PSFGauss(PSFModel): ...@@ -17,7 +17,7 @@ class PSFGauss(PSFModel):
self.fwhm = self.fwhmGauss(r=psfRa) self.fwhm = self.fwhmGauss(r=psfRa)
self.sigGauss = psfRa # 80% light radius self.sigGauss = psfRa # 80% light radius
self.sigSpin = sigSpin 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): def perfGauss(self, r, sig):
""" """
...@@ -122,4 +122,4 @@ class PSFGauss(PSFModel): ...@@ -122,4 +122,4 @@ class PSFGauss(PSFModel):
# return ell, beta, qr # return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians) 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
observation_sim/psf/PSFInterpSLS.py
View file @ 398f6272
...@@ -435,7 +435,16 @@ class PSFInterpSLS(PSFModel): ...@@ -435,7 +435,16 @@ class PSFInterpSLS(PSFModel):
# PSF_int_trans[ids_szero] = 0 # PSF_int_trans[ids_szero] = 0
# print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape) # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape)
PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans) 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 # from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans) # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
......
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