SLS bug fixed; saturation update; field distortion update

Catalog/C3Catalog.py

@@ -199,6 +199,8 @@ class C3Catalog(CatalogBase):
             stars = star_cat[str(pix)]
             self._load_gals(gals, pix_id=pix)
             self._load_stars(stars, pix_id=pix)
+            del gals
+            del stars
         print("number of objects in catalog: ", len(self.objs))
         del self.avGal
 
@@ -229,4 +231,6 @@ class C3Catalog(CatalogBase):
         # erg/s/cm2/A --> photo/s/m2/A
         all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
         sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        del wave
+        del flux
         return sed

ObservationSim/Instrument/Chip/Chip.py

@@ -320,6 +320,7 @@ class Chip(FocalPlane):
                 img.bounds, 
                 int(config["random_seeds"]["seed_flat"]))
             flat_normal = flat_img / np.mean(flat_img.array)
+            if self.survey_type == "photometric":
                 img *= flat_normal
             del flat_normal
             if config["output_setting"]["flat_output"] == False:
@@ -329,6 +330,7 @@ class Chip(FocalPlane):
         if config["ins_effects"]["shutter_effect"] == True:
             print("  Apply shutter effect", flush=True)
             shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3)    # shutter effect normalized image for this chip
+            if self.survey_type == "photometric":
                 img *= shuttimg
             if config["output_setting"]["shutter_output"] == True:    # output 16-bit shutter effect image with pixel value <=65535
                 shutt_gsimg = galsim.ImageUS(shuttimg*6E4)

ObservationSim/Instrument/Chip/Effects.py

@@ -3,7 +3,7 @@ from matplotlib.pyplot import flag
 import numpy as np
 from numpy.core.fromnumeric import mean, size
 from numpy.random import Generator, PCG64
-import math
+import math,copy
 from numba import jit
 from astropy import stats
 
@@ -87,10 +87,10 @@ def BadColumns(GSImage, seed=20240309, chipid=1):
     xposit = rgxpos.integers(low=int(xsize*0.05), high=int(xsize*0.95), size=(nbadsecA+nbadsecD))
     print(xposit+1)
     # signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1
-    if meanimg>0:
-        dn = rgdn.integers(low=meanimg*1.3+50, high=meanimg*2+150, size=(nbadsecA+nbadsecD)) #*signs
-    elif meanimg<0:
-        dn = rgdn.integers(low=meanimg*2-150, high=meanimg*1.3-50, size=(nbadsecA+nbadsecD)) #*signs
+    # if meanimg>0:
+    dn = rgdn.integers(low=np.abs(meanimg)*1.3+50, high=np.abs(meanimg)*2+150, size=(nbadsecA+nbadsecD)) #*signs
+    # elif meanimg<0:
+    #     dn = rgdn.integers(low=meanimg*2-150, high=meanimg*1.3-50, size=(nbadsecA+nbadsecD)) #*signs
     for badcoli in range(nbadsecA):
         GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli],1)))+dn[badcoli])
     for badcoli in range(nbadsecD):
@@ -106,6 +106,7 @@ def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=2021
         BiasLevel = np.zeros((nsecy,nsecx))
     elif bias_level>0:
         BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level
+    print(" Biases of 16 channels:\n",BiasLevel)
     arrshape = GSImage.array.shape
     secsize_x = int(arrshape[1]/nsecx)
     secsize_y = int(arrshape[0]/nsecy)
@@ -153,6 +154,22 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
     return GSImage
 
 
+def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
+    # Generate Gain non-uniformity, and multipy the different factors (mean~1 with sigma~1%) to the GS-Image
+    rg = Generator(PCG64(int(seed)))
+    Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1   # sigma~1%
+    Gain16 = Random16.reshape((nsecy,nsecx))/gain
+    print(seed-20210202, "Gains of 16 channels:\n", Gain16)
+    # arrshape = GSImage.array.shape
+    # secsize_x = int(arrshape[1]/nsecx)
+    # secsize_y = int(arrshape[0]/nsecy)
+    # for rowi in range(nsecy):
+    #     for coli in range(nsecx):
+    #         GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli]
+    # return GSImage
+    return Gain16
+
+
 def MakeFlatSmooth(GSBounds, seed):
     rg = Generator(PCG64(int(seed)))
     r1,r2,r3,r4 = rg.random(4)
@@ -235,104 +252,128 @@ def NonLinearity(GSImage, beta1=5E-7, beta2=0):
     return GSImage
 
 
-def chargeflow(ndarr, fullwell=10E4):
-    size_y,size_x = ndarr.shape
-    satpos_y = np.where(ndarr>=fullwell)[0]
-    satpos_x = np.where(ndarr>=fullwell)[1]
-    Nsatpix = len(satpos_y)
-
-    if Nsatpix==0:
-        # make no change for the image array
-        return ndarr
+########################################   Saturation & Bleeding Start    ###############################
 
+def BleedingTrail(aa, yy):
+    if aa<0.2:
+        aa=0.2
     else:
-        for i in range(Nsatpix):
-            ClumpFullwell=True
-            satcharge = ndarr[satpos_y[i],satpos_x[i]]-fullwell
-            ndarr[satpos_y[i],satpos_x[i]] = fullwell
-            
-            satpos_yi0 = satpos_y[i]
-
-            # Define the x,y=0,0 element of image array is the lower-left corner. So y decreases being downward.
-
-            # print('Charge Clump moves down')
-            chargedn = ((np.random.random()-0.5)*0.05+0.5)*satcharge
-            chargeup = satcharge - chargedn
-            fwi = 1
-            aa = np.log(chargedn/fullwell)**3*0.9   # blooming length begin to has e- less than fullwell
-            if aa < 0.05:
-                ClumpFullwell=False
-            # Test
+        pass
     try:
-                while ClumpFullwell==True:
-                    if satpos_y[i]<=0:
-                        break
-                    if ndarr[satpos_y[i]-1,satpos_x[i]]<fullwell:
-                        ndarr[satpos_y[i]-1,satpos_x[i]] = ndarr[satpos_y[i]-1,satpos_x[i]] + chargedn
-                        if ndarr[satpos_y[i]-1,satpos_x[i]]>=fullwell:
-                            fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa))
-                            if fx>5:
-                                fx=5
-                            faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa))
-                            rand_frac = 1-0.1*(fx-1)/(faa-1)
-
-                            chargedn = ndarr[satpos_y[i]-1,satpos_x[i]] - fullwell
-                            ndarr[satpos_y[i]-1,satpos_x[i]] = fullwell*rand_frac
-                            satpos_y[i] = satpos_y[i]-1
-                            if satpos_y[i]<0:
-                                ClumpFullwell=False
-                                break
-                        else:
-                            ClumpFullwell=False
-                        fwi += 1
-                    else:
-                        satpos_y[i] = satpos_y[i]-1
-                        if satpos_y[i]<0:
-                            ClumpFullwell=False
-                            break
+        fy = 0.5*(math.exp(math.log(yy+1)**3/aa)+np.exp(-1*math.log(yy+1)**3/aa))
+        faa= 0.5*(math.e+1/math.e)                
+        trail_frac = 1-0.1*(fy-1)/(faa-1)
+    except Exception as e:
+        print(e)
+        trail_frac = 1
+
+    return trail_frac
 
+def MakeTrail(imgarr, satuyxtuple, charge, fullwell=9e4, direction='up', trailcutfrac=0.9):
+    '''
+    direction: "up" or "down". For "up", bleeds along Y-decreasing direction; for "down", bleeds along Y-increasing direction.
+    '''
+    yi,xi = satuyxtuple
+    aa = np.log(charge/fullwell)**3              # scale length of the bleeding trail
+    yy = 1
 
-                # print('Charge Clump moves up')
-                ClumpFullwell=True
-                satpos_y[i] = satpos_yi0
-                fwi = 1
-                aa = np.log(chargeup/fullwell)**3*0.9   # blooming length at which it begins to have e- less than fullwell
-                if aa < 0.05:
-                    ClumpFullwell=False
-                while ClumpFullwell==True:
-                    if satpos_y[i]>=size_y-1:
+    while charge>0:
+        if yi<0 or yi>imgarr.shape[0]-1:
             break
-                    if ndarr[satpos_y[i]+1,satpos_x[i]]<fullwell:
-                        ndarr[satpos_y[i]+1,satpos_x[i]] = ndarr[satpos_y[i]+1,satpos_x[i]] + chargeup
-                        if ndarr[satpos_y[i]+1,satpos_x[i]]>=fullwell:
-                            fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa))
-                            if fx>5:
-                                fx=5
-                            faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa))
-                            rand_frac = 1-0.1*(fx-1)/(faa-1)
-                            chargeup = ndarr[satpos_y[i]+1,satpos_x[i]] - fullwell
-                            ndarr[satpos_y[i]+1,satpos_x[i]] = fullwell*rand_frac
-                            satpos_y[i] = satpos_y[i]+1
-                            if satpos_y[i]>=size_y-1:
-                                ClumpFullwell=False
+        if yi==0 or yi==imgarr.shape[0]-1:
+            imgarr[yi,xi] = fullwell
+            break
+        if direction=='up':
+            if imgarr[yi-1,xi]>=fullwell:
+                imgarr[yi,xi] = fullwell
+                yi-=1
+                continue
+        elif direction=='down':
+            if imgarr[yi+1,xi]>=fullwell:
+                imgarr[yi,xi] = fullwell
+                yi+=1
+                continue
+        if aa<=1:
+            while imgarr[yi,xi] >= fullwell:
+                imgarr[yi,xi] = fullwell
+                if direction=='up':
+                    imgarr[yi-1,xi] += charge
+                    charge = imgarr[yi-1,xi]-fullwell
+                    yi-=1
+                    if yi<0:
+                        break
+                elif direction=='down':
+                    imgarr[yi+1,xi] += charge
+                    charge = imgarr[yi+1,xi]-fullwell
+                    yi+=1
+                    if yi>imgarr.shape[0]:
                         break
         else:
-                            ClumpFullwell=False
-                        fwi += 1
+            # calculate bleeding trail:
+            trail_frac = BleedingTrail(aa,yy)
+
+            # put charge upwards
+            if trail_frac>=0.99:
+                imgarr[yi,xi] = fullwell
+                if direction=='up':
+                    yi-=1
+                elif direction=='down':
+                    yi+=1
+                yy += 1
             else:
-                        satpos_y[i] = satpos_y[i]+1
-                        if satpos_y[i]>=size_y-1:
-                            ClumpFullwell=False
+                if trail_frac<trailcutfrac:
                     break
-            except Exception as e:
-                print(e)
-                print(fwi, aa)
-                pass
+                charge = fullwell*trail_frac
+                imgarr[yi,xi] += charge
+                if imgarr[yi,xi]>fullwell:
+                    imgarr[yi,xi] = fullwell
+
+                if direction=='up':
+                    yi-=1
+                elif direction=='down':
+                    yi+=1
+                yy += 1
+
+    return imgarr
+
+
+def ChargeFlow(imgarr, fullwell=9E4):
+    size_y,size_x = imgarr.shape
+    satupos_y,satupos_x = np.where(imgarr>fullwell)
 
-    return ndarr
+    if satupos_y.shape[0]==0:
+        # make no change for the image array
+        return imgarr
+    elif satupos_y.shape[0]/imgarr.size > 0.5:
+        imgarr.fill(fullwell)
+        return imgarr
+
+    chargedict = {}
+    imgarrorig = copy.deepcopy(imgarr)
+
+    for yi,xi in zip(satupos_y,satupos_x):
+        yxidx = ''.join([str(yi),str(xi)])
+        chargedict[yxidx] = imgarrorig[yi,xi]-fullwell
 
+    for yi,xi in zip(satupos_y,satupos_x):
+        yxidx = ''.join([str(yi),str(xi)])
+        satcharge = chargedict[yxidx]
+        chargeup = ((np.random.random()-0.5)*0.05+0.5)*satcharge
+        chargedn = satcharge - chargeup
 
-def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=10e4):
+        try:
+            # Charge Clump moves up
+            if yi>=0 and yi<imgarr.shape[0]:
+                imgarr = MakeTrail(imgarr, (yi,xi), chargeup, fullwell=9e4, direction='up', trailcutfrac=0.9)
+                # Charge Clump moves down
+                imgarr = MakeTrail(imgarr, (yi,xi), chargedn, fullwell=9e4, direction='down', trailcutfrac=0.9)
+        except Exception as e:
+            print(e,'@pix ',(yi+1,xi+1))
+            return imgarr
+        
+    return imgarr
+
+def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=9e4):
     """
     To simulate digital detector's saturation and blooming effect. The blooming is along the read-out direction, perpendicular to the charge transfer direction. Charge clumpy overflows the pixel well will flow to two oposite directions with nearly same charges.
     Work together with chargeflow() function.
@@ -351,12 +392,14 @@ def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=10e4):
         for j in range(nsect_x):
             subimg = imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)]
 
-            subimg = chargeflow(subimg, fullwell=fullwell)
+            subimg = ChargeFlow(subimg, fullwell=fullwell)
 
             imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)] = subimg
 
     return GSImage
 
+#################################      Saturation & Bleeding End    ####################################
+
 
 def readout16(GSImage, rowi=0, coli=0, overscan_value=0):
     # readout image as 16 outputs of sub-images plus prescan & overscan.
@@ -672,7 +715,7 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2
     return CRmap.astype(np.int32), cr_event_size
 
 
-def ShutterEffectArr(GSImage, t_shutter=1.3, dist_bearing=735, dt=1E-3):
+def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E-3):
     # Generate Shutter-Effect normalized image
     # t_shutter: time of shutter movement
     # dist_bearing: distance between two bearings of shutter leaves
@@ -709,7 +752,10 @@ def ShutterEffectArr(GSImage, t_shutter=1.3, dist_bearing=735, dt=1E-3):
         s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb))
         brt[(idx>s1idx[i]) & (idx<s2idx[i])] += dt
 
-    brt = brt*2+(150-t_shutter*2)
+    if t_exp>t_shutter*2:
+        brt = brt*2+(t_exp-t_shutter*2)
+    else:
+        brt = brt*2
 
     x = (x-dist_bearing/2)*100
 

ObservationSim/MockObject/SkybackgroundMap.py

@@ -13,9 +13,12 @@ import os
 
 ###calculate sky map by sky SED
 
-def calculateSkyMap_split_g(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500, skyfn='param/skybackground/sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0,
+def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='param/skybackground/sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0,
                             split_pos=3685):
-    skyMap = np.ones([yLen, xLen], dtype='float32')
+    # skyMap = np.ones([yLen, xLen], dtype='float32')
+    #
+    # if isAlongY == 1:
+    #     skyMap = np.ones([xLen, yLen], dtype='float32')
 
     # for i in range(len(conf)):
     #     conf[i] = os.path.join(SLSSIM_PATH, conf[i])
@@ -24,9 +27,6 @@ def calculateSkyMap_split_g(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
     if np.size(conf) == 2:
         conf2 = conf[1]
 
-    if isAlongY == 1:
-        skyMap = np.ones([xLen, yLen], dtype='float32')
-
     skyImg = galsim.Image(skyMap, xmin=0, ymin=0)
 
     tbstart = blueLimit
@@ -68,6 +68,8 @@ def calculateSkyMap_split_g(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
             ssImg = galsim.ImageF(img_s)
             ssImg.setOrigin(origin_order_x, origin_order_y)
             bounds = ssImg.bounds & fImg.bounds
+            if bounds.area() == 0:
+                continue
             fImg[bounds] = fImg[bounds] + ssImg[bounds]
         
 
@@ -100,6 +102,8 @@ def calculateSkyMap_split_g(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
             ssImg = galsim.ImageF(img_s)
             ssImg.setOrigin(origin_order_x, origin_order_y)
             bounds = ssImg.bounds & fImg.bounds
+            if bounds.area() == 0:
+                continue
             fImg[bounds] = fImg[bounds] + ssImg[bounds]
 
         

ObservationSim/ObservationSim.py

@@ -10,6 +10,7 @@ from datetime import datetime
 from ObservationSim.Config import config_dir, ChipOutput
 from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
 from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
+from ObservationSim.Instrument.Chip import Effects
 from ObservationSim.MockObject import calculateSkyMap_split_g
 from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
 from ObservationSim._util import get_shear_field, makeSubDir_PointingList
@@ -89,10 +90,26 @@ class Observation(object):
         chip.img.wcs = wcs_fp
         if chip.survey_type == "photometric":
             sky_map = None
+        # elif chip.survey_type == "spectroscopic":
+        #     sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0)
         elif chip.survey_type == "spectroscopic":
-            sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0)
+            flat_normal = np.ones_like(chip.img.array)
+            if self.config["ins_effects"]["flat_fielding"] == True:
+                print("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID, flush=True)
+                print(chip.img.bounds, flush=True)
+                flat_img = Effects.MakeFlatSmooth(
+                    chip.img.bounds,
+                    int(self.config["random_seeds"]["seed_flat"]))
+                flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array)
+            if self.config["ins_effects"]["shutter_effect"] == True:
+                print("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID, flush=True)
+                shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735,
+                                                    dt=1E-3)  # shutter effect normalized image for this chip
+                flat_normal = flat_normal*shuttimg
+                flat_normal = np.array(flat_normal,dtype='float32')
+            sky_map = calculateSkyMap_split_g(skyMap=flat_normal, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0)
+            del flat_normal
 
-        # if pointing_type == 'MS':
         if pointing.pointing_type == 'MS':
             # Load catalogues and templates
             self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir)

ObservationSim/PSF/FieldDistortion.py

@@ -17,9 +17,12 @@ class FieldDistortion(object):
             self.fdModel = fdModel
 
     def isContainObj_FD(self, chip, pos_img):
-        ifdModel = self.fdModel["ccd" + chip.getChipLabel(chipID=chip.chipID)]["wave1"]
+        # ifdModel = self.fdModel["ccd" + chip.getChipLabel(chipID=chip.chipID)]["wave1"]
+        # x, y = pos_img.x, pos_img.y
+        # xLowI, xUpI, yLowI, yUpI = ifdModel["InterpLim"]
+        ifdModel = self.fdModel["wave1"]
+        xLowI, xUpI, yLowI, yUpI = ifdModel["interpLimit"]
         x, y = pos_img.x, pos_img.y
-        xLowI, xUpI, yLowI, yUpI = ifdModel["InterpLim"]
         if (xLowI - x)*(xUpI - x) > 0 or (yLowI - y) * (yUpI - y) > 0:
             return False
         return True
@@ -40,8 +43,18 @@ class FieldDistortion(object):
         """
         if not self.isContainObj_FD(chip=chip, pos_img=pos_img):
             return galsim.PositionD(-1, -1)
-        ifdModel = self.fdModel["ccd" + chip.getChipLabel(chipID=chip.chipID)]["wave1"]
+        # ifdModel = self.fdModel["ccd" + chip.getChipLabel(chipID=chip.chipID)]["wave1"]
+        ifdModel = self.fdModel["wave1"]
+        irsModel = self.fdModel["wave1"]["residual"]["ccd" + chip.getChipLabel(chipID=chip.chipID)]
         x, y = pos_img.x, pos_img.y
         x = ifdModel["xImagePos"](x, y)[0][0]
         y = ifdModel["yImagePos"](x, y)[0][0]
+        x1LowI, x1UpI, y1LowI, y1UpI = irsModel["interpLimit"]
+        if (x1LowI-x)*(x1UpI-x) <=0 and (y1LowI-y)*(y1UpI-y)<=0:
+            dx = irsModel["xResidual"](x, y)[0][0]
+            dy = irsModel["yResidual"](x, y)[0][0]
+            x += dx
+            y += dy
+        else:
+            return galsim.PositionD(-1, -1)
         return galsim.PositionD(x, y)

run_C3.pbs

@@ -7,13 +7,13 @@
 ##mpdboot -n 10 -f  ./mpd.hosts
 
 ##PBS -j oe
-#PBS -l nodes=comput112:ppn=80
+#PBS -l nodes=comput108:ppn=80
 #####PBS -q longq
 #PBS -q batch
 #PBS -u fangyuedong
 
-NP=80
+NP=40
 date
 echo $NP
 
-mpirun -np $NP python /public/home/fangyuedong/20211203/CSST/run_sim.py config_C3.yaml -c /public/home/fangyuedong/20211203/CSST/config
+mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_C3.yaml -c /public/home/fangyuedong/test/CSST/config