diff --git a/ObservationSim/Config/Config.py b/ObservationSim/Config/Config.py
index 74837df4fcc5e42ff506a39a9dae2ab572c4d53b..3d6724b139c2fc2f6cad907185d4947e04d34329 100755
--- a/ObservationSim/Config/Config.py
+++ b/ObservationSim/Config/Config.py
@@ -20,6 +20,7 @@ def config_dir(config, work_dir=None, data_dir=None):
# PSF data directory
if config["psf_setting"]["psf_model"] == "Interp":
path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"])
+ path_dict["psf_sls_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_sls_dir"])
return path_dict
diff --git a/ObservationSim/Instrument/Chip/ChipUtils.py b/ObservationSim/Instrument/Chip/ChipUtils.py
index 0a7b804e1097265e94cadc7479ca0808f73fd9f7..381d310da490cf46704635b74f737b7d755074f4 100644
--- a/ObservationSim/Instrument/Chip/ChipUtils.py
+++ b/ObservationSim/Instrument/Chip/ChipUtils.py
@@ -21,6 +21,22 @@ def log_info(msg, logger=None):
else:
print(msg, flush=True)
+def getChipSLSGratingID(chipID):
+ gratingID = ['','']
+ if chipID == 1: gratingID = ['GI2', 'GI1']
+ if chipID == 2: gratingID = ['GV4', 'GV3']
+ if chipID == 3: gratingID = ['GU2', 'GU1']
+ if chipID == 4: gratingID = ['GU4', 'GU3']
+ if chipID == 5: gratingID = ['GV2', 'GV1']
+ if chipID == 10: gratingID = ['GI4', 'GI3']
+ if chipID == 21: gratingID = ['GI6', 'GI5']
+ if chipID == 26: gratingID = ['GV8', 'GV7']
+ if chipID == 27: gratingID = ['GU6', 'GU5']
+ if chipID == 28: gratingID = ['GU8', 'GU7']
+ if chipID == 29: gratingID = ['GV6', 'GV5']
+ if chipID == 30: gratingID = ['GI8', 'GI7']
+ return gratingID
+
def getChipSLSConf(chipID):
confFile = ''
if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py
index 81ce517cb7bc4fbcefea1905dac16afaf63e2322..81c29cd63ff9e0109194e072bba9cd24d6fe00c2 100755
--- a/ObservationSim/MockObject/Galaxy.py
+++ b/ObservationSim/MockObject/Galaxy.py
@@ -248,10 +248,27 @@ class Galaxy(MockObject):
xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
+
+ branges = np.zeros([len(bandpass_list), 2])
+
+ # print(hasattr(psf_model, 'bandranges'))
+
+ if hasattr(psf_model, 'bandranges'):
+ if psf_model.bandranges is None:
+ return 2, None
+ if len(psf_model.bandranges) != len(bandpass_list):
+ return 2, None
+ branges = psf_model.bandranges
+ else:
+ for i in range(len(bandpass_list)):
+ branges[i, 0] = bandpass_list[i].blue_limit * 10
+ branges[i, 1] = bandpass_list[i].red_limit * 10
+
for i in range(len(bandpass_list)):
- bandpass = bandpass_list[i]
+ # bandpass = bandpass_list[i]
+ brange = branges[i]
- psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
+ # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
@@ -272,14 +289,14 @@ class Galaxy(MockObject):
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
- gal = galsim.Convolve(psf, gal)
+ # gal = galsim.Convolve(psf, gal)
- if not big_galaxy: # Not apply PSF for very big galaxy
- gal = galsim.Convolve(psf, gal)
- # if fd_shear is not None:
- # gal = gal.shear(fd_shear)
+ # if not big_galaxy: # Not apply PSF for very big galaxy
+ # gal = galsim.Convolve(psf, gal)
+ # # if fd_shear is not None:
+ # # gal = gal.shear(fd_shear)
- starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset)
+ starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
@@ -301,12 +318,16 @@ class Galaxy(MockObject):
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
+ psf_model=psf_model, bandNo=i + 1,
+ grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
@@ -318,12 +339,16 @@ class Galaxy(MockObject):
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
+ psf_model=psf_model, bandNo=i + 1,
+ grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp_p1
del sdp_p2
@@ -331,25 +356,33 @@ class Galaxy(MockObject):
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
+ psf_model=psf_model, bandNo=i + 1,
+ grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
elif grating_split_pos_chip>=gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
+ psf_model=psf_model, bandNo=i + 1,
+ grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin)
- del psf
+ # del psf
return 1, pos_shear
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py
index ce6e3ba62a9148350610ed64b7faeab19611a79f..8e361b9b2a79250f2bc05fc0c162382be581be6f 100755
--- a/ObservationSim/MockObject/MockObject.py
+++ b/ObservationSim/MockObject/MockObject.py
@@ -4,7 +4,7 @@ import astropy.constants as cons
from astropy import wcs
from astropy.table import Table
-from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders
+from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders, convolveImg
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, \
getABMAG
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
@@ -223,9 +223,69 @@ class MockObject(object):
del stamp
del spec_orders
+ def addSLStoChipImageWithPSF(self, sdp=None, chip=None, pos_img_local = [1,1], psf_model=None, bandNo = 1, grating_split_pos=3685, local_wcs=None, pos_img=None):
+ spec_orders = sdp.compute_spec_orders()
+ for k, v in spec_orders.items():
+ img_s = v[0]
+ # print(bandNo,k)
+ try:
+ psf, pos_shear = psf_model.get_PSF(chip, pos_img_local = pos_img_local, bandNo = bandNo, galsimGSObject=True, g_order = k, grating_split_pos=grating_split_pos)
+ except:
+ psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
+
+ psf_img = psf.drawImage(nx=100, ny=100, wcs = local_wcs)
+
+ psf_img_m = psf_img.array
+
+ #########################################################
+ # DEBUG
+ #########################################################
+ # ids_p = psf_img_m < 0
+ # psf_img_m[ids_p] = 0
+
+ # from astropy.io import fits
+ # fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m)
+
+ # print("DEBUG: orig_off is", orig_off)
+ nan_ids = np.isnan(img_s)
+ if img_s[nan_ids].shape[0] > 0:
+ img_s[nan_ids] = 0
+ print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
+ #########################################################
+ img_s, orig_off = convolveImg(img_s, psf_img_m)
+ origin_order_x = v[1] - orig_off[0]
+ origin_order_y = v[2] - orig_off[1]
+
+
+ specImg = galsim.ImageF(img_s)
+ # photons = galsim.PhotonArray.makeFromImage(specImg)
+ # photons.x += origin_order_x
+ # photons.y += origin_order_y
+
+ # xlen_imf = int(specImg.xmax - specImg.xmin + 1)
+ # ylen_imf = int(specImg.ymax - specImg.ymin + 1)
+ # stamp = galsim.ImageF(xlen_imf, ylen_imf)
+ # stamp.wcs = local_wcs
+ # stamp.setOrigin(origin_order_x, origin_order_y)
+
+ specImg.wcs = local_wcs
+ specImg.setOrigin(origin_order_x, origin_order_y)
+
+ bounds = specImg.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
+ if bounds.area() == 0:
+ continue
+ chip.img.setOrigin(0, 0)
+ chip.img[bounds] = chip.img[bounds] + specImg[bounds]
+ # stamp[bounds] = chip.img[bounds]
+ # # chip.sensor.accumulate(photons, stamp)
+ # chip.img[bounds] = stamp[bounds]
+ chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
+ # del stamp
+ del spec_orders
+ return pos_shear
+
def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
-
if normFilter is not None:
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
@@ -263,22 +323,38 @@ class MockObject(object):
xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062, 'C': 4.035447379743442,
'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
+
+ branges = np.zeros([len(bandpass_list),2])
+
+ if hasattr(psf_model,'bandranges'):
+ if psf_model.bandranges is None:
+ return 2, None
+ if len(psf_model.bandranges) != len(bandpass_list):
+ return 2, None
+ branges = psf_model.bandranges
+ else:
+ for i in range(len(bandpass_list)):
+ branges[i, 0] = bandpass_list[i].blue_limit * 10
+ branges[i, 1] = bandpass_list[i].red_limit * 10
+
for i in range(len(bandpass_list)):
- bandpass = bandpass_list[i]
- psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
- folding_threshold=folding_threshold)
+ # bandpass = bandpass_list[i]
+ brange = branges[i]
+
+ # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
+ # folding_threshold=folding_threshold)
star = galsim.DeltaFunction(gsparams=gsp)
star = star.withFlux(tel.pupil_area * exptime)
- star = galsim.Convolve(psf, star)
+ psf_tmp = galsim.Gaussian(sigma=0.002)
+ star = galsim.Convolve(psf_tmp, star)
- starImg = star.drawImage(nx=100, ny=100, wcs=chip_wcs_local, offset=offset)
+ starImg = star.drawImage(nx=60, ny=60, wcs=chip_wcs_local, offset=offset)
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0,0)
gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
-
if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse
@@ -293,12 +369,15 @@ class MockObject(object):
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local = [xcenter_p1,ycenter_p1],
+ psf_model=psf_model, bandNo = i+1, grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
@@ -310,12 +389,15 @@ class MockObject(object):
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
+ psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp_p1
del sdp_p2
@@ -323,23 +405,29 @@ class MockObject(object):
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
+ psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
- band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
+ band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
- self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+ pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
+ psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
+ local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
- del psf
+ # del psf
return 1, pos_shear
def SNRestimate(self, img_obj, flux, noise_level=0.0, seed=31415):
diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py
index 3f3f52b01a1d9b061099491d5e49f677979411ee..7c8b51efa9f39b07e485f00506a53eaff163f118 100644
--- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py
+++ b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py
@@ -155,7 +155,7 @@ class SpecDisperser(object):
sensitivity_beam = ysens
len_spec_x = len(dx)
- len_spec_y = int(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0]) + 1)
+ len_spec_y = int(abs(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0])) + 1)
beam_sh = (self.img_sh[0] + len_spec_y, self.img_sh[1] + len_spec_x)
modelf = zeros(product(beam_sh), dtype=float)
diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py
index 6f51570bd9426fd40c94b019f17914cde6f65150..02aee7700f411d53bf05850bfee514d3f4898f45 100755
--- a/ObservationSim/MockObject/_util.py
+++ b/ObservationSim/MockObject/_util.py
@@ -571,4 +571,14 @@ def convolveGaussXorders(img=None, sigma = 1):
convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
return convImg, offset
+def convolveImg(img=None, psf = None):
+ from astropy.modeling.models import Gaussian2D
+ from scipy import signal
+
+ convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
+ offset_x = int(psf.shape[1]/2. + 0.5) - 1
+ offset_y = int(psf.shape[0]/2. + 0.5) - 1
+ offset = [offset_x,offset_y]
+ return convImg, offset
+
diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py
index e35d700886658b13b415f2910d95bbd17f2c0489..9500e473bfa57738c2902be20e335acd0ce7e33f 100755
--- a/ObservationSim/ObservationSim.py
+++ b/ObservationSim/ObservationSim.py
@@ -15,7 +15,7 @@ from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensio
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from ObservationSim.Instrument.Chip import Effects
from ObservationSim.Straylight import calculateSkyMap_split_g
-from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
+from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS
from ObservationSim._util import get_shear_field, makeSubDir_PointingList
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
@@ -63,7 +63,10 @@ class Observation(object):
if self.config["psf_setting"]["psf_model"] == "Gauss":
psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"])
elif self.config["psf_setting"]["psf_model"] == "Interp":
- psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"])
+ if chip.survey_type == "spectroscopic":
+ psf_model = PSFInterpSLS(chip, filt,PSF_data_prefix=self.path_dict["psf_sls_dir"])
+ else:
+ psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"])
else:
chip_output.Log_error("unrecognized PSF model type!!", flush=True)
@@ -198,6 +201,10 @@ class Observation(object):
obj = self.cat.objs[j]
+ # (DEBUG)
+ # if obj.getMagFilter(filt)>20:
+ # continue
+
# load and convert SED; also caculate object's magnitude in all CSST bands
try:
sed_data = self.cat.load_sed(obj)
diff --git a/ObservationSim/PSF/PSFInterpSLS.py b/ObservationSim/PSF/PSFInterpSLS.py
new file mode 100644
index 0000000000000000000000000000000000000000..7dc1624c2ae4305cac78650044e46c95e6de57a2
--- /dev/null
+++ b/ObservationSim/PSF/PSFInterpSLS.py
@@ -0,0 +1,570 @@
+'''
+PSF interpolation for CSST-Sim
+
+NOTE: [iccd, iwave, ipsf] are counted from 1 to n, but [tccd, twave, tpsf] are counted from 0 to n-1
+'''
+
+import sys
+import time
+import copy
+import numpy as np
+import scipy.spatial as spatial
+import galsim
+import h5py
+
+from ObservationSim.PSF.PSFModel import PSFModel
+from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
+import os
+from astropy.io import fits
+
+from astropy.modeling.models import Gaussian2D
+from scipy import signal
+
+
+LOG_DEBUG = False #***#
+NPSF = 900 #***# 30*30
+PixSizeInMicrons = 5. #***# in microns
+
+
+###find neighbors-KDtree###
+# def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
+# """
+# find nearest neighbors by 2D-KDTree
+#
+# Parameters:
+# tx, ty (float, float): a given position
+# px, py (numpy.array, numpy.array): position data for tree
+# dr (float-optional): distance
+# dn (int-optional): nearest-N
+# OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
+#
+# Returns:
+# dataq (numpy.array): index
+# """
+# datax = px
+# datay = py
+# tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
+#
+# dataq=[]
+# rr = dr
+# if OnlyDistance == True:
+# dataq = tree.query_ball_point([tx, ty], rr)
+# if OnlyDistance == False:
+# while len(dataq) < dn:
+# dataq = tree.query_ball_point([tx, ty], rr)
+# rr += dr
+# dd = np.hypot(datax[dataq]-tx, datay[dataq]-ty)
+# ddSortindx = np.argsort(dd)
+# dataq = np.array(dataq)[ddSortindx[0:dn]]
+# return dataq
+#
+# ###find neighbors-hoclist###
+# def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
+# if np.max(partx) > nhocx*dhocx:
+# print('ERROR')
+# sys.exit()
+# if np.max(party) > nhocy*dhocy:
+# print('ERROR')
+# sys.exit()
+#
+# npart = partx.size
+# hoclist= np.zeros(npart, dtype=np.int32)-1
+# hoc = np.zeros([nhocy, nhocx], dtype=np.int32)-1
+# for ipart in range(npart):
+# ix = int(partx[ipart]/dhocx)
+# iy = int(party[ipart]/dhocy)
+# hoclist[ipart] = hoc[iy, ix]
+# hoc[iy, ix] = ipart
+# return hoc, hoclist
+#
+# def hocFind(px, py, dhocx, dhocy, hoc, hoclist):
+# ix = int(px/dhocx)
+# iy = int(py/dhocy)
+#
+# neigh=[]
+# it = hoc[iy, ix]
+# while it != -1:
+# neigh.append(it)
+# it = hoclist[it]
+# return neigh
+#
+# def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None):
+# nhocy = nhocx = 20
+#
+# pxMin = np.min(px)
+# pxMax = np.max(px)
+# pyMin = np.min(py)
+# pyMax = np.max(py)
+#
+# dhocx = (pxMax - pxMin)/(nhocx-1)
+# dhocy = (pyMax - pyMin)/(nhocy-1)
+# partx = px - pxMin +dhocx/2
+# party = py - pyMin +dhocy/2
+#
+# if hoc is None:
+# hoc, hoclist = hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy)
+# return hoc, hoclist
+#
+# if hoc is not None:
+# tx = tx - pxMin +dhocx/2
+# ty = ty - pyMin +dhocy/2
+# itx = int(tx/dhocx)
+# ity = int(ty/dhocy)
+#
+# ps = [-1, 0, 1]
+# neigh=[]
+# for ii in range(3):
+# for jj in range(3):
+# ix = itx + ps[ii]
+# iy = ity + ps[jj]
+# if ix < 0:
+# continue
+# if iy < 0:
+# continue
+# if ix > nhocx-1:
+# continue
+# if iy > nhocy-1:
+# continue
+#
+# #neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist)
+# it = hoc[iy, ix]
+# while it != -1:
+# neigh.append(it)
+# it = hoclist[it]
+# #neigh.append(neightt)
+# #ll = [i for k in neigh for i in k]
+# if dn != -1:
+# ptx = np.array(partx[neigh])
+# pty = np.array(party[neigh])
+# dd = np.hypot(ptx-tx, pty-ty)
+# idx = np.argsort(dd)
+# neigh= np.array(neigh)[idx[0:dn]]
+# return neigh
+#
+#
+# ###PSF-IDW###
+# def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False):
+# """
+# psf interpolation by IDW
+#
+# Parameters:
+# px, py (float, float): position of the target
+# PSFMat (numpy.array): image
+# cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers
+# IDWindex (int-optional): the power index of IDW
+# OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation
+#
+# Returns:
+# psfMaker (numpy.array)
+# """
+#
+# minimum_psf_weight = 1e-8
+# ref_col = px
+# ref_row = py
+#
+# ngy, ngx = PSFMat[0, :, :].shape
+# npsf = PSFMat[:, :, :].shape[0]
+# psfWeight = np.zeros([npsf])
+#
+# if OnlyNeighbors == True:
+# if hoc is None:
+# neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False)
+# if hoc is not None:
+# neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist)
+#
+# neighFlag = np.zeros(npsf)
+# neighFlag[neigh] = 1
+#
+# for ipsf in range(npsf):
+# if OnlyNeighbors == True:
+# if neighFlag[ipsf] != 1:
+# continue
+#
+# dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
+# if IDWindex == 1:
+# psfWeight[ipsf] = dist
+# if IDWindex == 2:
+# psfWeight[ipsf] = dist**2
+# if IDWindex == 3:
+# psfWeight[ipsf] = dist**3
+# if IDWindex == 4:
+# psfWeight[ipsf] = dist**4
+# psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
+# psfWeight[ipsf] = 1./psfWeight[ipsf]
+# psfWeight /= np.sum(psfWeight)
+#
+# psfMaker = np.zeros([ngy, ngx], dtype=np.float32)
+# for ipsf in range(npsf):
+# if OnlyNeighbors == True:
+# if neighFlag[ipsf] != 1:
+# continue
+#
+# iPSFMat = PSFMat[ipsf, :, :].copy()
+# ipsfWeight = psfWeight[ipsf]
+#
+# psfMaker += iPSFMat * ipsfWeight
+# psfMaker /= np.nansum(psfMaker)
+#
+# return psfMaker
+
+
+
+###define PSFInterp###
+class PSFInterpSLS(PSFModel):
+ def __init__(self, chip, filt,PSF_data_prefix="", sigSpin=0, psfRa=0.15, pix_size = 0.005):
+ if LOG_DEBUG:
+ print('===================================================')
+ print('DEBUG: psf module for csstSim ' \
+ +time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
+ print('===================================================')
+
+ self.sigSpin = sigSpin
+ self.sigGauss = psfRa
+ self.grating_ids = chip_utils.getChipSLSGratingID(chip.chipID)
+ _,self.grating_type = chip.getChipFilter(chipID=chip.chipID)
+ self.data_folder = PSF_data_prefix
+ self.getPSFDataFromFile(filt)
+ self.pixsize = pix_size # um
+
+ def getPSFDataFromFile(self, filt):
+ gratingInwavelist = {'GU':0,'GV':1,'GI':2}
+ grating_orders = ['0','1']
+ waveListFn = self.data_folder + '/wavelist.dat'
+ wavelists = np.loadtxt(waveListFn)
+ self.waveList = wavelists[:,gratingInwavelist[self.grating_type]]
+ bandranges = np.zeros([4,2])
+ midBand = (self.waveList[0:3] + self.waveList[1:4])/2.*10000.
+ bandranges[0,0] = filt.blue_limit
+ bandranges[1:4,0] = midBand
+ bandranges[0:3, 1] = midBand
+ bandranges[3,1] = filt.red_limit
+
+ self.bandranges = bandranges
+
+ self.grating1_data = {}
+ g_folder = self.data_folder + '/' + self.grating_ids[0] + '/'
+ for g_order in grating_orders:
+ g_folder_order = g_folder + 'PSF_Order_' + g_order + '/'
+ grating_order_data = {}
+ for bandi in [1,2,3,4]:
+ subBand_data = {}
+ subBand_data['bandrange'] = bandranges[bandi-1]
+ final_folder = g_folder_order + str(bandi) + '/'
+ print(final_folder)
+ pca_fs = os.listdir(final_folder)
+ for fname in pca_fs:
+ if ('_PCs.fits' in fname) and (fname[0] != '.'):
+ fname_ = final_folder + fname
+ hdu = fits.open(fname_)
+ subBand_data['band_data'] = hdu
+ grating_order_data['band'+str(bandi)] = subBand_data
+ self.grating1_data['order'+g_order] = grating_order_data
+
+ self.grating2_data = {}
+ g_folder = self.data_folder + '/' + self.grating_ids[1] + '/'
+ for g_order in grating_orders:
+ g_folder_order = g_folder + 'PSF_Order_' + g_order + '/'
+ grating_order_data = {}
+ for bandi in [1, 2, 3, 4]:
+ subBand_data = {}
+ subBand_data['bandrange'] = bandranges[bandi - 1]
+ final_folder = g_folder_order + str(bandi) + '/'
+ print(final_folder)
+ pca_fs = os.listdir(final_folder)
+ for fname in pca_fs:
+ if ('_PCs.fits' in fname) and (fname[0] != '.'):
+ fname_ = final_folder + fname
+ hdu = fits.open(fname_)
+ subBand_data['band_data'] = hdu
+ grating_order_data['band' + str(bandi)] = subBand_data
+ self.grating2_data['order' + g_order] = grating_order_data
+
+ #
+ #
+ #
+ # def _getPSFwave(self, iccd, PSF_data_file, PSF_data_prefix):
+ # # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
+ # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
+ # nwave = len(fq.keys())
+ # fq.close()
+ # return nwave
+ #
+ #
+ # def _loadPSF(self, iccd, PSF_data_file, PSF_data_prefix):
+ # psfSet = []
+ # # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
+ # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
+ # for ii in range(self.nwave):
+ # iwave = ii+1
+ # psfWave = []
+ #
+ # fq_iwave = fq['w_{:}'.format(iwave)]
+ # for jj in range(self.npsf):
+ # ipsf = jj+1
+ # psfInfo = {}
+ # psfInfo['wavelength']= fq_iwave['wavelength'][()]
+ #
+ # fq_iwave_ipsf = fq_iwave['psf_{:}'.format(ipsf)]
+ # psfInfo['pixsize'] = PixSizeInMicrons
+ # psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()]
+ # psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()]
+ # psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()]
+ # psfInfo['image_y'] = fq_iwave_ipsf['image_y'][()]
+ # psfInfo['centroid_x']= fq_iwave_ipsf['cx'][()]
+ # psfInfo['centroid_y']= fq_iwave_ipsf['cy'][()]
+ # psfInfo['psfMat'] = fq_iwave_ipsf['psfMat'][()]
+ #
+ # psfWave.append(psfInfo)
+ # psfSet.append(psfWave)
+ # fq.close()
+ #
+ # if LOG_DEBUG:
+ # print('psfSet has been loaded:', flush=True)
+ # print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True)
+ # return psfSet
+ #
+ #
+ # def _findWave(self, bandpass):
+ # if isinstance(bandpass,int):
+ # twave = bandpass
+ # return twave
+ #
+ # for twave in range(self.nwave):
+ # bandwave = self.PSF_data[twave][0]['wavelength']
+ # if bandpass.blue_limit < bandwave and bandwave < bandpass.red_limit:
+ # return twave
+ # return -1
+ #
+ #
+
+ def convolveWithGauss(self, img=None, sigma=1):
+
+ offset = int(np.ceil(sigma * 3))
+ g_size = 2 * offset + 1
+ m_cen = int(g_size / 2)
+ print('-----',g_size)
+ g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma)
+ yp, xp = np.mgrid[0:g_size, 0:g_size]
+ g_PSF = g_PSF_(xp, yp)
+ psf = g_PSF / g_PSF.sum()
+ convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
+ convImg = convImg/np.sum(convImg)
+ return convImg
+
+
+ def get_PSF(self, chip, pos_img_local = [1000,1000], bandNo = 1, galsimGSObject=True, folding_threshold=5.e-3, g_order = 'A', grating_split_pos=3685):
+ """
+ Get the PSF at a given image position
+
+ Parameters:
+ chip: A 'Chip' object representing the chip we want to extract PSF from.
+ pos_img: A 'galsim.Position' object representing the image position.
+ bandpass: A 'galsim.Bandpass' object representing the wavelength range.
+ pixSize: The pixels size of psf matrix
+ findNeighMode: 'treeFind' or 'hoclistFind'
+ Returns:
+ PSF: A 'galsim.GSObject'.
+ """
+ order_IDs = {'A': '1', 'B': '0' ,'C': '0', 'D': '0', 'E': '0'}
+ contam_order_sigma = {'C':0.28032344707964174,'D':0.39900182912061344,'E':1.1988309797685412} #arcsec
+ x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
+ y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
+ # print(pos_img.x - x_start)
+ pos_img_x = pos_img_local[0] + x_start
+ pos_img_y = pos_img_local[1] + y_start
+ pos_img = galsim.PositionD(pos_img_x, pos_img_y)
+ if pos_img_local[0] < grating_split_pos:
+ psf_data = self.grating1_data
+ else:
+ psf_data = self.grating2_data
+
+ grating_order = order_IDs[g_order]
+ # if grating_order in ['-2','-1','2']:
+ # grating_order = '1'
+
+
+ # if grating_order in ['0', '1']:
+ psf_order = psf_data['order'+grating_order]
+ psf_order_b = psf_order['band'+str(bandNo)]
+ psf_b_dat = psf_order_b['band_data']
+ pos_p = psf_b_dat[1].data
+ pc_coeff = psf_b_dat[2].data
+ pcs = psf_b_dat[0].data
+ # print(max(pos_p[:,0]), min(pos_p[:,0]),max(pos_p[:,1]), min(pos_p[:,1]))
+ # print(chip.x_cen, chip.y_cen)
+ # print(pos_p)
+ px = pos_img.x*chip.pix_size
+ py = pos_img.y*chip.pix_size
+
+ dist2=(pos_p[:,1] - px)*(pos_p[:,1] - px) + (pos_p[:,0] - py)*(pos_p[:,0] - py)
+ temp_sort_dist = np.zeros([dist2.shape[0],2])
+ temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0],1)
+ temp_sort_dist[:, 1] = dist2
+ # print(temp_sort_dist)
+ dits2_sortlist = sorted(temp_sort_dist, key=lambda x:x[1])
+ # print(dits2_sortlist)
+ nearest4p = np.zeros([4,2])
+ pc_coeff_4p = np.zeros([pc_coeff.data.shape[0],4])
+
+ for i in np.arange(4):
+ smaller_ids = int(dits2_sortlist[i][0])
+ nearest4p[i, 0] = pos_p[smaller_ids, 1]
+ nearest4p[i, 1] = pos_p[smaller_ids, 0]
+ pc_coeff_4p[:,i] = pc_coeff[:,smaller_ids]
+ idw_dist = 1/(np.sqrt((px-nearest4p[:,0]) * (px-nearest4p[:,0]) + (py-nearest4p[:,1]) * (py-nearest4p[:,1])))
+
+ coeff_int = np.zeros(pc_coeff.data.shape[0])
+ for i in np.arange(4):
+ coeff_int = coeff_int + pc_coeff_4p[:,i]*idw_dist[i]
+ coeff_int = coeff_int / np.sum(coeff_int)
+
+ npc = 10
+ m_size = int(pcs.shape[0]**0.5)
+ PSF_int = np.dot(pcs[:,0:npc],coeff_int[0:npc]).reshape(m_size,m_size)
+
+ # PSF_int = PSF_int/np.sum(PSF_int)
+ PSF_int_trans = np.flipud(np.fliplr(PSF_int))
+ PSF_int_trans = np.fliplr(PSF_int_trans.T)
+ # PSF_int_trans = np.abs(PSF_int_trans)
+ # ids_szero = PSF_int_trans<0
+ # PSF_int_trans[ids_szero] = 0
+ # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape)
+ PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans)
+
+ # from astropy.io import fits
+ # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
+
+
+ # if g_order in ['C','D','E']:
+ # g_simgma = contam_order_sigma[g_order]/pixel_size_arc
+ # PSF_int_trans = self.convolveWithGauss(PSF_int_trans,g_simgma)
+ # n_m_size = int(m_size/2)
+ #
+ # n_PSF_int = np.zeros([n_m_size, n_m_size])
+ #
+ # for i in np.arange(n_m_size):
+ # for j in np.arange(n_m_size):
+ # n_PSF_int[i,j] = np.sum(PSF_int[2*i:2*i+2, 2*j:2*j+2])
+ #
+ # n_PSF_int = n_PSF_int/np.sum(n_PSF_int)
+
+ # chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
+ # chip.img.wcs = galsim.wcs.AffineTransform
+ if galsimGSObject:
+
+ # imPSFt = np.zeros([257,257])
+ # imPSFt[0:256, 0:256] = imPSF
+ # # imPSFt[120:130, 0:256] = 1.
+
+ pixel_size_arc = np.rad2deg(self.pixsize * 1e-3 / 28) * 3600
+ img = galsim.ImageF(PSF_int_trans, scale=pixel_size_arc)
+ gsp = galsim.GSParams(folding_threshold=folding_threshold)
+ ############TEST: START
+ # Use sheared PSF to test the PSF orientation
+ # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
+ ############TEST: END
+ self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
+ # if g_order in ['C','D','E']:
+ # add_psf = galsim.Gaussian(sigma=contam_order_sigma[g_order], flux=1.0)
+ # self.psf = galsim.Convolve(self.psf, add_psf)
+ wcs = chip.img.wcs.local(pos_img)
+ scale = galsim.PixelScale(0.074)
+ self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
+
+ # return self.PSFspin(x=px/0.01, y=py/0.01)
+ return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
+
+ return PSF_int_trans, PSF_int
+
+
+ # pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize
+ #
+ # # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID'
+ # twave = self._findWave(bandpass)
+ # if twave == -1:
+ # print("!!!PSF bandpass does not match.")
+ # exit()
+ # PSFMat = self.psfMat[twave]
+ # cen_col= self.cen_col[twave]
+ # cen_row= self.cen_row[twave]
+ #
+ # px = (pos_img.x - chip.cen_pix_x)*0.01
+ # py = (pos_img.y - chip.cen_pix_y)*0.01
+ # if findNeighMode == 'treeFind':
+ # imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, PSFCentroidWgt=True)
+ # if findNeighMode == 'hoclistFind':
+ # assert(self.hoc != 0), 'hoclist should be built correctly!'
+ # imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)
+ #
+ # ############TEST: START
+ # TestGaussian = False
+ # if TestGaussian:
+ # gsx = galsim.Gaussian(sigma=0.04)
+ # #pointing_pa = -23.433333
+ # imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array
+ # ############TEST: END
+ #
+ # if galsimGSObject:
+ # imPSFt = np.zeros([257,257])
+ # imPSFt[0:256, 0:256] = imPSF
+ # # imPSFt[120:130, 0:256] = 1.
+ #
+ # img = galsim.ImageF(imPSFt, scale=pixSize)
+ # gsp = galsim.GSParams(folding_threshold=folding_threshold)
+ # ############TEST: START
+ # # Use sheared PSF to test the PSF orientation
+ # # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
+ # ############TEST: END
+ # self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
+ # wcs = chip.img.wcs.local(pos_img)
+ # scale = galsim.PixelScale(0.074)
+ # self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
+ #
+ # # return self.PSFspin(x=px/0.01, y=py/0.01)
+ # return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
+ # return imPSF
+ #
+ # def PSFspin(self, x, y):
+ # """
+ # The PSF profile at a given image position relative to the axis center
+ #
+ # Parameters:
+ # theta : spin angles in a given exposure in unit of [arcsecond]
+ # dx, dy: relative position to the axis center in unit of [pixels]
+ #
+ # Return:
+ # Spinned PSF: g1, g2 and axis ratio 'a/b'
+ # """
+ # a2Rad = np.pi/(60.0*60.0*180.0)
+ #
+ # ff = self.sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
+ # rc = np.sqrt(x*x + y*y)
+ # cpix = rc*(self.sigSpin*a2Rad)
+ #
+ # beta = (np.arctan2(y,x) + np.pi/2)
+ # ell = cpix**2/(2.0*ff**2+cpix**2)
+ # qr = np.sqrt((1.0+ell)/(1.0-ell))
+ # PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
+ # return self.psf.shear(PSFshear), PSFshear
+
+from ObservationSim.Instrument import Filter, FilterParam, Chip
+import yaml
+if __name__ == '__main__':
+ configfn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/config/config_C6_dev.yaml'
+ with open(configfn, "r") as stream:
+ try:
+ config = yaml.safe_load(stream)
+ for key, value in config.items():
+ print (key + " : " + str(value))
+ except yaml.YAMLError as exc:
+ print(exc)
+ chip = Chip(chipID=1,config=config)
+ filter_id, filter_type = chip.getChipFilter()
+ filt = Filter(filter_id=filter_id,
+ filter_type=filter_type,
+ filter_param=FilterParam())
+
+ psf_i = PSFInterpSLS(chip, filt,PSF_data_prefix="/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/")
+ pos_img = galsim.PositionD(x=25155, y=-22060)
+ psf_im = psf_i.get_PSF(chip, pos_img = pos_img, g_order = '1')
+
diff --git a/ObservationSim/PSF/__init__.py b/ObservationSim/PSF/__init__.py
index d0ea11850ebb2728337a2b2ce76d714b4591fb5a..80bf7d3685fc1aba1b88aa8a4535163ca39d55fa 100755
--- a/ObservationSim/PSF/__init__.py
+++ b/ObservationSim/PSF/__init__.py
@@ -2,4 +2,5 @@ from .PSFModel import PSFModel
from .PSFGauss import PSFGauss
# from .PSFInterp.PSFInterp import PSFInterp
from .PSFInterp import PSFInterp
+from .PSFInterpSLS import PSFInterpSLS
from .FieldDistortion import FieldDistortion
\ No newline at end of file
diff --git a/config/config_C6_dev.yaml b/config/config_C6_dev.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..0619cb6cb80297200eecec248e94b792fe8649a7
--- /dev/null
+++ b/config/config_C6_dev.yaml
@@ -0,0 +1,228 @@
+---
+###############################################
+#
+# 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
+<<<<<<< HEAD
+work_dir: "/share/home/zhangxin/CSST_SIM/CSST_new_sim/csst-simulation/"
+=======
+work_dir: "/share/home/fangyuedong/new_sim/workplace/"
+# work_dir: "/share/C6_new_sim_2sq"
+>>>>>>> new_sim
+data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
+run_name: "C6_new_sim_2sq_run1"
+project_cycle: 6
+run_counter: 1
+
+# 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: YES
+
+ # 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: "/share/simudata/CSSOSDataProductsSims/data/"
+ 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: [0]
+
+ # 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: True
+
+ # 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
+ mag_sat_margin: -15.
+
+ # 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: "/share/simudata/CSSOSDataProductsSims/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: YES # 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: YES # 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: YES # Whether to add bad columns
+ add_hotpixels: YES # Whether to add hot pixels
+ add_deadpixels: YES # 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
+...