bug fixes & extrapolate PSF for saturated objects

config/obs_config_SCI.yaml

@@ -62,10 +62,10 @@ call_sequence:
     hot_pixels: YES
     dead_pixels: YES
     bad_columns: YES
-  # Apply response nonlinearity
-  nonlinearity: {}
   # Apply CCD Saturation & Blooming
   blooming: {}
+  # Apply response nonlinearity
+  nonlinearity: {}
   # Run CTE simulation
   # CTE_effect: {}
   # Add prescan and overscan

observation_sim/instruments/chip/libBF/diffusion_X1.c

@@ -37,7 +37,7 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi
   {
     printf("Adding BF effect...\n");
     //setup BF correlation fliter
-    float neX;
+    float neX, neXtemp;
     float neP1 = 50000;
     float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302};
     float neP2 = 10000;
@@ -56,15 +56,18 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi
         neX = arr_ima[j+i*ny];
         if(neX >= 10000)
         {
+          neXtemp = neX;
+          if(neXtemp > 100000)
+              neXtemp = 100000;
           bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0;
-          bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575;
-          bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652;
-          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335;
-          bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]);
-          bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118;
-          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]);
-          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083;
-          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043;
+          bfa[0][1]=bfa[0][-1]=linearInterp(neXtemp, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575;
+          bfa[-1][0]=bfa[1][0]=linearInterp(neXtemp, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652;
+          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neXtemp, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335;
+          bfa[0][-2]=bfa[0][2]=linearInterp(neXtemp, neP1, bfaP1[4], neP2, bfaP2[4]);
+          bfa[-2][0]=bfa[2][0]=linearInterp(neXtemp, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118;
+          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neXtemp, neP1, bfaP1[6], neP2, bfaP2[6]);
+          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neXtemp, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083;
+          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neXtemp, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043;
         }
         else
         {

observation_sim/mock_objects/Galaxy.py

@@ -115,7 +115,7 @@ class Galaxy(MockObject):
 
         # Set Galsim Parameters
         if self.getMagFilter(filt) <= 15 and (not big_galaxy):
-            folding_threshold = 5.e-4
+            folding_threshold = 5.e-8
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
@@ -170,8 +170,11 @@ class Galaxy(MockObject):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             # Get PSF model
+            EXTRA = False
+            if self.getMagFilter(filt) <= filt.mag_saturation-1.:
+                EXTRA = True
             psf, pos_shear = psf_model.get_PSF(
-                chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
+                chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold, extrapolate=EXTRA)
 
             if self.bfrac == 0:
                 gal_temp = disk

observation_sim/mock_objects/MockObject.py

@@ -122,7 +122,7 @@ class MockObject(object):
             return 2, None
         # Set Galsim Parameters
         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)
@@ -160,14 +160,17 @@ class MockObject(object):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             # Get PSF model
+            EXTRA = False
+            if self.getMagFilter(filt) <= filt.mag_saturation-1.:
+                EXTRA = True
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
-                                               folding_threshold=folding_threshold)
+                                               folding_threshold=folding_threshold, extrapolate=EXTRA)
             # star = galsim.DeltaFunction(gsparams=gsp)
             # star = star.withFlux(nphotons)
             # star = galsim.Convolve(psf, star)
             star = psf.withFlux(nphotons)
 
-            stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
+            stamp = star.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset)
             if np.sum(np.isnan(stamp.array)) > 0:
                 continue
             stamp.setCenter(x_nominal, y_nominal)

observation_sim/psf/PSFInterp.py

@@ -13,6 +13,7 @@ import galsim
 import h5py
 
 from observation_sim.psf.PSFModel import PSFModel
+from observation_sim.psf._util import psf_extrapolate
 
 
 NPSF = 900  # ***# 30*30
@@ -324,7 +325,7 @@ class PSFInterp(PSFModel):
                 return twave
         return -1
 
-    def get_PSF(self, chip, pos_img, bandpass, galsimGSObject=True, findNeighMode='treeFind', folding_threshold=5.e-3, pointing_pa=0.0):
+    def get_PSF(self, chip, pos_img, bandpass, galsimGSObject=True, findNeighMode='treeFind', folding_threshold=5.e-3, pointing_pa=0.0, extrapolate=False, ngg=2048):
         """
         Get the PSF at a given image position
 
@@ -358,20 +359,17 @@ class PSFInterp(PSFModel):
             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 extrapolate is True:
+            ccdList = [ 6,  7,  8,  9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23, 24, 25]
+            rr_trim_list = [72, 64, 96, 88, 64, 72, 72, 76, 72, 72, 76, 72, 72, 64, 88, 96, 64, 72]
+            imPSF = psf_extrapolate(imPSF, rr_trim=rr_trim_list[ccdList.index(chip.chipID)], ngg=ngg)
 
         if galsimGSObject:
-            imPSFt = np.zeros([257, 257])
-            imPSFt[0:256, 0:256] = imPSF
-            # imPSFt[120:130, 0:256] = 1.
+            if extrapolate is True:
+                imPSFt = imPSF
+            else:
+                imPSFt = np.zeros([257, 257])
+                imPSFt[0:256, 0:256] = imPSF
 
             img = galsim.ImageF(imPSFt, scale=pixSize)
             gsp = galsim.GSParams(folding_threshold=folding_threshold)

observation_sim/psf/_util.py

+import numpy as np
+from scipy.interpolate import interp1d
+
+def binningPSF(img, ngg):
+    imgX = img.reshape(ngg, img.shape[0]//ngg, ngg, img.shape[1]//ngg).mean(-1).mean(1)
+    return imgX
+
+
+def radial_average_at_pixel(image, center_x, center_y, dr=10):
+    # Get coordinates relative to the specified center pixel (x, y)
+    y, x = np.indices(image.shape)
+    r = np.sqrt((x - center_x)**2 + (y - center_y)**2)
+
+    # Set up bins
+    max_radius = int(r.max())  # Maximum distance from the center pixel
+    radial_bins = np.arange(0, max_radius, dr)
+
+    # Compute average value in each bin
+    radial_means = []
+    for i in range(len(radial_bins) - 1):
+        mask = (r >= radial_bins[i]) & (r < radial_bins[i + 1])
+        if np.any(mask):
+            radial_means.append(image[mask].mean())
+        else:
+            radial_means.append(0)  # In case no pixels are in the bin
+    return radial_bins[:-1], radial_means  # Exclude last bin since no mean calculated
+
+
+def psf_extrapolate(psf, rr_trim=64, ngg=256):
+    # ngg = 256
+    # extrapolate PSF
+    if True:
+        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.log(radii[1:])
+        means_log = np.log(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_temp + np.finfo(float).eps)
+        arr[rim  > rr_trim] = 0
+        arr[arr==0] = y_new[arr==0]
+        psf = np.exp(arr)
+        psf[rim  > int(ngg/2)] = 0
+    imPSF = psf  # binningPSF(psf, int(ngg/2))
+    imPSF = imPSF/np.nansum(imPSF)
+    return imPSF