Merge branch 'new_sim' of https://csst-tb.bao.ac.cn/code/csst_sim/csst-simulation into new_sim

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
 

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']

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)):
-            bandpass = bandpass_list[i]
+                branges[i, 0] = bandpass_list[i].blue_limit * 10
+                branges[i, 1] = bandpass_list[i].red_limit * 10
 
-            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
+        for i in range(len(bandpass_list)):
+            # 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)
             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.):

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)):
-            bandpass = bandpass_list[i]
-            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
-                                               folding_threshold=folding_threshold)
+                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]
+            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):

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)

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
+
 

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,6 +63,9 @@ 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":
+            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)

ObservationSim/PSF/PSFInterpSLS.py

+'''
+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

@@ -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

config/config_C6_dev.yaml

+---
+###############################################
+#
+#  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
+...