add led flat to image and sls; add brighter source psf for sls

catalog/C9_Catalog.py

@@ -470,8 +470,9 @@ class Catalog(CatalogBase):
             self.objs.append(obj)
 
     def free_mem(self, **kward):
-        self.starDDL.freeGlobeData()
-        del self.starDDL
+        if not self.config["catalog_options"]["galaxy_only"]:
+            self.starDDL.freeGlobeData()
+            del self.starDDL
 
     def _load(self, **kwargs):
         self.objs = []

observation_sim/ObservationSim.py

@@ -62,6 +62,8 @@ class Observation(object):
         # Get flat, shutter, and PRNU images
         chip.flat_img, _ = chip_utils.get_flat(
             img=chip.img, seed=int(self.config["random_seeds"]["seed_flat"]))
+        if chip.chipID <= 30:
+            chip.flat_img = chip.flat_img*chip_utils.get_innerflat(chip = chip)
         if chip.chipID > 30:
             chip.shutter_img = np.ones_like(chip.img.array)
         else:

observation_sim/instruments/chip/chip_utils.py

@@ -172,6 +172,12 @@ def get_flat(img, seed):
     flat_normal = flat_img / np.mean(flat_img.array)
     return flat_img, flat_normal
 
+def get_innerflat(chip = None, filt = None):
+    from observation_sim.mock_objects import FlatLED
+    led_obj = FlatLED(chip, filt)
+    flat_img = led_obj.getInnerFlat()
+    return flat_img
+
 
 def add_cosmic_rays(img, chip, exptime=150, seed=0):
     cr_map, cr_event_num = effects.produceCR_Map(

observation_sim/mock_objects/FlatLED.py

@@ -50,6 +50,8 @@ fluxLED = {'LED1': 15, 'LED2': 15, 'LED3': 12.5, 'LED4': 9, 'LED5': 9,
 #            'LED14': 10}
 
 mirro_eff = {'GU': 0.61, 'GV': 0.8, 'GI': 0.8}
+
+bandtoLed = {'NUV':['LED1','LED2'], 'u':['LED13','LED14'], 'g':['LED3','LED4','LED5'], 'r':['LED6','LED7'], 'i':['LED8'], 'z':['LED9','LED10'], 'y':['LED10'], 'GU':['LED1','LED2','LED13','LED14'], 'GV':['LED3','LED4','LED5','LED6'], 'GI':['LED7','LED8','LED9','LED10']}
 # mirro_eff = {'GU':1, 'GV':1, 'GI':1}
 
 
@@ -69,10 +71,19 @@ class FlatLED(MockObject):
                 with pkg_resources.path('observation_sim.mock_objects.data.led', "") as ledDir:
                     self.flatDir = ledDir.as_posix()
 
+
+    def getInnerFlat(self):
+        ledflats = bandtoLed[self.chip.filter_type]
+        iFlat = np.zeros([self.chip.npix_y, self.chip.npix_x])
+        for nled in ledflats:
+            iFlat = iFlat + self.getLEDImage(led_type=nled, LED_Img_flag =False)
+        iFlat = iFlat / len(ledflats)
+        return iFlat
+
     ###
     # return LED flat, e/s
     ###
-    def getLEDImage(self, led_type='LED1'):
+    def getLEDImage(self, led_type='LED1', LED_Img_flag =True):
         # cwave = cwaves[led_type]
         flat = fits.open(os.path.join(self.flatDir, 'model_' +
                          cwaves_name[led_type] + 'nm.fits'))
@@ -102,7 +113,10 @@ class FlatLED(MockObject):
             N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]),
             method='linear')
         U = U/np.mean(U)
-        flatImage = U*fluxLED[led_type]*1000
+
+        flatImage = U
+        if LED_Img_flag:
+            flatImage = flatImage*fluxLED[led_type]*1000
         gc.collect()
         return flatImage
 

observation_sim/mock_objects/Galaxy.py

@@ -335,18 +335,21 @@ class Galaxy(MockObject):
                 pos_img_local[1] = pos_img.y - y_start
                 nnx = 0
                 nny = 0
-                for order in ["A", "B"]:
+                for order in ["B", "A"]:
+                    EXTRA = False
+                    if self.getMagFilter(filt) <= filt.mag_saturation-2.:
+                        EXTRA = True
                     psf, pos_shear = psf_model.get_PSF(
-                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
+                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos, extrapolate = EXTRA, ngg=3072)
                     star_p = galsim.Convolve(psf, gal)
                     if nnx == 0:
                         galImg = star_p.drawImage(
-                            wcs=chip_wcs_local, offset=offset)
+                            wcs=chip_wcs_local, offset=offset, method='no_pixel')
                         nnx = galImg.xmax - galImg.xmin + 1
                         nny = galImg.ymax - galImg.ymin + 1
                     else:
                         galImg = star_p.drawImage(
-                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset)
+                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset, method='no_pixel')
                     galImg.setOrigin(0, 0)
                     # n1 = np.sum(np.isinf(galImg.array))
                     # n2 = np.sum(np.isnan(galImg.array))

observation_sim/mock_objects/MockObject.py

@@ -386,7 +386,7 @@ class MockObject(object):
             sedNormFactor = 1.
 
         if self.getMagFilter(filt) <= 15:
-            folding_threshold = 5.e-4
+            folding_threshold = 5.e-8
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
@@ -445,19 +445,24 @@ class MockObject(object):
                 pos_img_local[1] = pos_img.y - y_start
                 nnx = 0
                 nny = 0
-                for order in ["A", "B"]:
+                for order in ["B", "A"]:
+                    EXTRA = False
+                    if self.getMagFilter(filt) <= filt.mag_saturation-2.:
+                        EXTRA = True
+                        nnx = 2048
+                        nny = 2048
                     psf, pos_shear = psf_model.get_PSF(
-                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
+                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos, extrapolate=EXTRA, ngg=3072)
                     # star_p = galsim.Convolve(psf, star)
                     star_p = psf.withFlux(tel.pupil_area * exptime)
                     if nnx == 0:
                         starImg = star_p.drawImage(
-                            wcs=chip_wcs_local, offset=offset)
+                            wcs=chip_wcs_local, offset=offset, method='no_pixel')
                         nnx = starImg.xmax - starImg.xmin + 1
                         nny = starImg.ymax - starImg.ymin + 1
                     else:
                         starImg = star_p.drawImage(
-                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset)
+                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset, method='no_pixel')
                     # n1 = np.sum(np.isinf(starImg.array))
                     # n2 = np.sum(np.isnan(starImg.array))
                     # if n1>0 or n2 > 0:
@@ -470,7 +475,7 @@ class MockObject(object):
                 psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
                 # star_p = galsim.Convolve(psf, star)
                 star_p = psf.withFlux(tel.pupil_area * exptime)
-                starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
+                starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset, method='no_pixel')
                 starImg.setOrigin(0, 0)
                 for order in ["A", "B", "C", "D", "E"]:
                     starImg_List.append(starImg)

observation_sim/psf/PSFInterpSLS.py

@@ -23,6 +23,9 @@ from astropy.modeling.models import Gaussian2D
 from scipy import signal, interpolate
 import datetime
 import gc
+from astropy.io import fits
+
+from observation_sim.psf._util import psf_extrapolate, psf_extrapolate1
 # from jax import numpy as jnp
 
 LOG_DEBUG = False  # ***#
@@ -354,7 +357,7 @@ class PSFInterpSLS(PSFModel):
         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):
+    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, extrapolate=False, ngg=2048):
         """
         Get the PSF at a given image position
 
@@ -435,6 +438,9 @@ class PSFInterpSLS(PSFModel):
         # PSF_int_trans[ids_szero] = 0
         # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape)
         PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans)
+        PSF_int_trans = PSF_int_trans-np.min(PSF_int_trans)
+        PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans)
+        # fits.writeto('/home/zhangxin/CSST_SIM/CSST_sim_develop/psf_test/psf.fits',PSF_int_trans)
         # DEBGU
         ids_szero = PSF_int_trans < 0
         n01 = PSF_int_trans[ids_szero].shape[0]
@@ -444,7 +450,16 @@ class PSFInterpSLS(PSFModel):
         if n1 > 0 or n2 > 0:
             print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d" %
                   (n1, n2, n01))
-
+        if extrapolate is True:
+            # for rep_i in np.arange(0, 2, 1):
+                # PSF_int_trans[rep_i,:] = 1e9*pow(10,rep_i)
+                # PSF_int_trans[-1-rep_i,:]  = 1e9*pow(10,rep_i)
+                # PSF_int_trans[:,rep_i] = 1e9*pow(10,rep_i)
+                # PSF_int_trans[:,-1-rep_i] = 1e9*pow(10,rep_i)
+
+            PSF_int_trans = psf_extrapolate1(PSF_int_trans, ngg=ngg)
+            # fits.writeto('/home/zhangxin/CSST_SIM/CSST_sim_develop/psf_test/psf_large.fits',PSF_int_trans) 
+        
         ####
         # from astropy.io import fits
         # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)

observation_sim/psf/_util.py

@@ -63,3 +63,58 @@ def psf_extrapolate(psf, rr_trim=64, ngg=256):
     imPSF = psf  # binningPSF(psf, int(ngg/2))
     imPSF = imPSF/np.nansum(imPSF)
     return imPSF
+
+def psf_extrapolate1(psf, rr_trim=64, ngg=256):
+    # ngg = 256
+    # extrapolate PSF
+    if True:
+        psf_enlar = np.zeros([ngg, ngg])
+        psf_enlar[int(ngg/2-128):int(ngg/2+128), int(ngg/2-128):int(ngg/2+128)] = psf
+        xim = np.arange(ngg)-ngg/2
+        xim, yim = np.meshgrid(xim, xim)
+        rim = np.sqrt(xim**2 + yim**2)
+
+        psf_temp = psf_enlar
+        # psf_temp[rim >= rr_trim] = 0
+        psf_temp[rim >= ngg/2-2] = np.finfo(float).eps
+        radii, means = radial_average_at_pixel(psf_temp, ngg/2, ngg/2, dr=2)
+
+        radii_log = np.log(radii[1:])
+        # radii_log = radii[1:]
+        means_log = np.log(means[1:])
+
+
+        # xim = np.arange(256)-128
+        # xim, yim = np.meshgrid(xim, xim)
+        # rim = np.sqrt(xim**2 + yim**2)
+
+        # # rr_trim  = 96
+        # psf_temp = psf
+        # psf_temp[rim > rr_trim] = 0
+        # radii, means = radial_average_at_pixel(psf_temp, 128, 128, dr=4)
+
+        # radii_log = np.log10(radii[1:])
+        # # radii_log = radii[1:]
+        # means_log = np.log10(means[1:])
+
+        finite_mask = np.isfinite(means_log)
+        f_interp = interp1d(radii_log[finite_mask][:-1], means_log[finite_mask][:-1], kind='linear', fill_value="extrapolate")
+
+        # ngg = 1024
+        # xim = np.arange(ngg)-int(ngg/2)
+        # xim, yim = np.meshgrid(xim, xim)
+        # rim = np.sqrt(xim**2 + yim**2)
+        # rim[int(ngg/2), int(ngg/2)] = np.finfo(float).eps  # 1e-7
+        rim_log = np.log(rim)
+        y_new = f_interp(rim_log)
+
+        arr = np.zeros([ngg, ngg])
+        arr[int(ngg/2-128):int(ngg/2+128), int(ngg/2-128):int(ngg/2+128)] = np.log(psf + np.finfo(float).eps)
+        arr[rim > 128-2] = 0
+        arr[arr == 0] = y_new[arr == 0]
+        psf_n = np.exp(arr)
+        # psf_n[int(ngg/2-128):int(ngg/2+128), int(ngg/2-128):int(ngg/2+128)] = psf
+        # psf[rim > int(ngg/2)] = 0
+    imPSF = psf_n  # binningPSF(psf, int(ngg/2))
+    # imPSF = imPSF/np.nansum(imPSF)
+    return imPSF