fixed photometry/bad column issue, modify detector effects order

ObservationSim/Config/ChipOutput.py

@@ -37,8 +37,10 @@ class ChipOutput(object):
         fmt5 = "%10s %8.4f %8.4f %8.4f\n"
         self.hdr = hdr1 + hdr2 + hdr3 + hdr4 + hdr5
         self.fmt = fmt1 + fmt2 + fmt3 + fmt4 + fmt5
+        print("pointing_type = %s\n"%(pointing_type))
         if pointing_type == 'MS':
             self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
+            print("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
             self.cat.write(self.hdr)
 
     def updateHDR(self, hdr):

ObservationSim/Instrument/Chip/Chip.py

@@ -30,16 +30,13 @@ class Chip(FocalPlane):
         self.read_noise = 5.0   # e/pix
         self.dark_noise = 0.02  # e/pix/s
         self.pix_scale  = 0.074 # pixel scale
-        # self.gain       = 1.0
-        # self.bias_level = 1000  # e-/pix
         self.gain = float(config["gain"])
         self.bias_level = float(config["bias_level"])
         self.overscan   = 1000
-        self.exptime    = 150   # second
-        # self.dark_exptime = 300  # second
-        # self.flat_exptime = 150  # second
+        self.exptime    = float(config["exp_time"])   # second
         self.dark_exptime = float(config['dark_exptime'])
         self.flat_exptime = float(config['flat_exptime'])
+        self.readout_time = float(config['readout_time'])
 
         # A chip ID must be assigned
         self.chipID = int(chipID)
@@ -234,60 +231,6 @@ class Chip(FocalPlane):
         noise = self.dark_noise * exptime + self.read_noise**2
         return noise
 
-
-    def addNoise_phot(self, img, exptime=150.0, sky_noise=0., seed=31415):
-        rng = galsim.BaseDeviate(seed)
-        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
-        img.addNoise(dark_noise)
-        ccd_noise = galsim.CCDNoise(rng, sky_level=sky_noise, gain=self.gain, read_noise=self.read_noise)
-        img.addNoise(ccd_noise)
-        return img
-
-    def addNoise_spec(self, config, tel, img, sky_map, exptime=150.0, seed=31415):
-        if img.array.shape != sky_map.shape:
-            raise ValueError("The shape img and sky_map must be equal.")
-
-        # n_img_arrar = (img.array + sky_map) * tel.pupil_area * exptime
-        # Should be the following?
-        rng = galsim.BaseDeviate(seed)
-        noise_img = galsim.Image(img, copy=True)
-        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
-        img.addNoise(dark_noise)
-        if config["abs_back"].lower() == "y":
-            img += (sky_map * tel.pupil_area * exptime)
-        ccd_noise = galsim.CCDNoise(rng, gain=self.gain, read_noise=self.read_noise)
-        img.addNoise(ccd_noise)
-        return img
-        # pNoise = self.dark_noise * exptime
-        # noise_img = galsim.Image(n_img_arrar, copy=True)
-        # dNose = galsim.PoissonNoise(rng=rng, sky_level=pNoise)
-        # noise_img.addNoise(dNose)
-        # rdNoise = galsim.GaussianNoise(rng=rng, sigma=self.read_noise)
-        # noise_img.addNoise(rdNoise)
-        # return noise_img
-
-    def addNoise(self, config, tel, filt, img, sky_map=None, exptime=150.0):
-        if self.survey_type == "photometric":
-            sky_level = filt.getSkyNoise(exptime=exptime, gain=self.gain)
-            img = self.addNoise_phot(
-                img=img, 
-                exptime=exptime, 
-                sky_noise=filt.getSkyNoise(exptime=exptime, gain=self.gain), 
-                seed=int(config["seed_skynoise"]))
-            if config["abs_back"].lower() == "y":
-                img += sky_level
-        elif self.survey_type == "spectroscopic":
-            if sky_map is None:
-                raise ValueError("Must provide sky_map for spectroscopic chip")
-            img = self.addNoise_spec(
-                config=config, 
-                tel=tel, 
-                img=img, 
-                sky_map=sky_map, 
-                exptime=exptime, 
-                seed=int(config["seed_skynoise"]))
-        return img
-
     def getChipSLSConf(self):
         confFile = ''
         if self.chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
@@ -351,7 +294,7 @@ class Chip(FocalPlane):
         fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
         hdu1.writeto(fname, output_verify='ignore', overwrite=True)
 
-    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS'):
+    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None):
         SeedGainNonuni=int(config["seed_gainNonUniform"])
         SeedBiasNonuni=int(config["seed_biasNonUniform"])
         SeedRnNonuni = int(config["seed_rnNonUniform"])
@@ -368,6 +311,17 @@ class Chip(FocalPlane):
         else:
             BoolDeadPix = False
 
+        # Add sky background
+        if sky_map == None:
+            sky_map = filt.getSkyNoise(exptime=self.exptime)
+        elif img.array.shape != sky_map.shape:
+            raise ValueError("The shape img and sky_map must be equal.")
+        elif tel is not None: # If sky_map is given in flux
+            sky_map = sky_map * tel.pupil_area * self.exptime
+        if config["add_back"].lower() == "y":
+            img += sky_map
+        del sky_map
+
         # Apply flat-field large scale structure for one chip
         if config["flat_fielding"].lower() == "y":
             print("  Creating and applying Flat-Fielding", flush=True)
@@ -392,30 +346,19 @@ class Chip(FocalPlane):
                 del shutt_gsimg
             del shuttimg
 
-        # Apply PRNU effect and output PRNU flat file:
-        if config["prnu_effect"].lower() == "y":
-            print("  Applying PRNU effect", flush=True)
-            prnu_img = effects.PRNU_Img(
-                xsize=self.npix_x, 
-                ysize=self.npix_y, 
-                sigma=0.01, 
-                seed=int(config["seed_prnu"]+self.chipID))
-            img *= prnu_img
-            if config["prnu_output"].lower() == "y":
-                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
-            if config["flat_output"].lower() == "n":
-                del prnu_img
-
-        # Add dark current
-        img += self.dark_noise*exptime
+        # Add Poisson noise
+        seed = int(config["seed_poisson"]) + pointing_ID*30 + self.chipID
+        rng_poisson = galsim.BaseDeviate(seed)
+        poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.)
+        img.addNoise(poisson_noise)
 
         # Add cosmic-rays
         if config["cosmic_ray"].lower() == "y" and pointing_type=='MS':
             print("  Adding Cosmic-Ray", flush=True)
             cr_map = effects.produceCR_Map(
                 xLen=self.npix_x, yLen=self.npix_y, 
-                exTime=exptime, 
-                cr_pixelRatio=0.003, 
+                exTime=self.exptime+0.5*self.readout_time, 
+                cr_pixelRatio=0.003*(1+0.5*self.readout_time/self.exptime), 
                 gain=self.gain, 
                 attachedSizes=self.attachedSizes,
                 seed=SeedCosmicRay+pointing_ID*30+self.chipID)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
@@ -442,15 +385,78 @@ class Chip(FocalPlane):
                 exptime=150.)
             del crmap_gsimg
 
+        # Apply PRNU effect and output PRNU flat file:
+        if config["prnu_effect"].lower() == "y":
+            print("  Applying PRNU effect", flush=True)
+            prnu_img = effects.PRNU_Img(
+                xsize=self.npix_x, 
+                ysize=self.npix_y, 
+                sigma=0.01, 
+                seed=int(config["seed_prnu"]+self.chipID))
+            img *= prnu_img
+            if config["prnu_output"].lower() == "y":
+                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
+            if config["flat_output"].lower() == "n":
+                del prnu_img
+
+        # Add dark current
+        if config["add_dark"].lower() == "y":
+            dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(self.exptime+0.5*self.readout_time)))
+            img.addNoise(dark_noise)
+
+        # Add Hot Pixels or/and Dead Pixels
+        rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
+        badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
+        img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)
+
+        # Apply Bad lines 
+        if config["add_badcolumns"].lower() == "y":
+            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID)
+
         # Add Bias level
         if config["add_bias"].lower() == "y":
             print("  Adding Bias level and 16-channel non-uniformity")
-            # img += float(config["bias_level"])
             img = effects.AddBiasNonUniform16(img, 
                 bias_level=float(config["bias_level"]), 
                 nsecy = 2, nsecx=8, 
                 seed=SeedBiasNonuni+self.chipID)
 
+        # Apply Nonlinearity on the chip image
+        if config["non_linear"].lower() == "y":
+            print("  Applying Non-Linearity on the chip image", flush=True)
+            img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0)
+
+        # Apply CCD Saturation & Blooming
+        if config["saturbloom"].lower() == "y":
+            print("  Applying CCD Saturation & Blooming")
+            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)
+
+        # Apply CTE Effect
+        if config["cte_trail"].lower() == "y":
+            print("  Apply CTE Effect")
+            img = effects.CTE_Effect(GSImage=img, threshold=27)
+
+        # Add Read-out Noise
+        if config["add_readout"].lower() == "y":
+            seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID
+            rng_readout = galsim.BaseDeviate(seed)
+            readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
+            img.addNoise(readout_noise)
+
+
+        # Apply Gain & Quantization
+        print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
+        img = effects.ApplyGainNonUniform16(
+            img, gain=self.gain, 
+            nsecy = 2, nsecx=8, 
+            seed=SeedGainNonuni+self.chipID)
+        img.array[img.array > 65535] = 65535
+        img.replaceNegative(replace_value=0)
+        img.quantize()
+
+        ######################################################################################
+        # Output images for calibration pointing
+        ######################################################################################
         # Bias output
         if config["bias_output"].lower() == "y" and pointing_type=='CAL':
             print("  Output N frame Bias files", flush=True)
@@ -466,7 +472,7 @@ class Chip(FocalPlane):
                         cr_map = effects.produceCR_Map(
                             xLen=self.npix_x, yLen=self.npix_y, 
                             exTime=0.01, 
-                            cr_pixelRatio=0.003*0.01/self.exptime, 
+                            cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), 
                             gain=self.gain, 
                             attachedSizes=self.attachedSizes,
                             seed=SeedCosmicRay+pointing_ID*30+self.chipID+1)
@@ -477,7 +483,11 @@ class Chip(FocalPlane):
                 # Non-Linearity for Bias
                 if config["non_linear"].lower() == "y":
                     print("  Applying Non-Linearity on the Bias image", flush=True)
-                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=1.e-7, beta2=1.e-10)
+                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0)
+
+                # Apply Bad lines 
+                if config["add_badcolumns"].lower() == "y":
+                    BiasCombImg = effects.BadColumns(BiasCombImg-float(config["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID) + float(config["bias_level"])-5
 
                 BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
@@ -493,7 +503,7 @@ class Chip(FocalPlane):
                 datetime_obs = datetime.fromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
-                timestamp_obs += 5 * 60
+                timestamp_obs += 10 * 60
                 self.outputCal(
                     img=BiasCombImg,
                     ra_cen=ra_cen,
@@ -517,7 +527,7 @@ class Chip(FocalPlane):
             elif config["add_bias"].lower() == "n":
                 biaslevel = 0
                 overscan = 0
-            darklevel = self.dark_noise*self.flat_exptime
+            darklevel = self.dark_noise*(self.flat_exptime+0.5*self.readout_time)
             for i in range(NFlat):
                 FlatSingle = flat_img * prnu_img + darklevel
                 FlatCombImg, FlatTag = effects.MakeFlatNcomb(
@@ -526,15 +536,15 @@ class Chip(FocalPlane):
                     read_noise=self.read_noise,
                     gain=1, 
                     overscan=overscan, 
-                    biaslevel=biaslevel,
+                    biaslevel=0,
                     seed_bias=SeedDefective+self.chipID
                     )
                 if config["cosmic_ray"].lower() == "y":
                     if config["cray_differ"].lower() == "y":
                         cr_map = effects.produceCR_Map(
                             xLen=self.npix_x, yLen=self.npix_y, 
-                            exTime=self.flat_exptime, 
-                            cr_pixelRatio=0.003*self.flat_exptime/self.exptime, 
+                            exTime=self.flat_exptime+0.5*self.readout_time, 
+                            cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), 
                             gain=self.gain, 
                             attachedSizes=self.attachedSizes,
                             seed=SeedCosmicRay+pointing_ID*30+self.chipID+3)
@@ -544,19 +554,35 @@ class Chip(FocalPlane):
 
                 if config["non_linear"].lower() == "y":
                     print("  Applying Non-Linearity on the Flat image", flush=True)
-                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=1.e-7, beta2=1.e-10)
+                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0)
 
                 if config["cte_trail"].lower() == "y":
                     FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)
 
+                # Add Hot Pixels or/and Dead Pixels
+                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
+                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
+                FlatCombImg = effects.DefectivePixels(FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)
+
                 # Apply Bad lines 
                 if config["add_badcolumns"].lower() == "y":
                     FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID)
 
-                # Add Hot Pixels or/and Dead Pixels
-                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
-                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
-                FlatCombImg = effects.DefectivePixels(FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=self.bias_level)
+                # Add Bias level
+                if config["add_bias"].lower() == "y":
+                    print("  Adding Bias level and 16-channel non-uniformity")
+                    # img += float(config["bias_level"])
+                    FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, 
+                        bias_level=biaslevel, 
+                        nsecy = 2, nsecx=8, 
+                        seed=SeedBiasNonuni+self.chipID)
+                
+                # Add Read-out Noise
+                if config["add_readout"].lower() == "y":
+                    seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID
+                    rng_readout = galsim.BaseDeviate(seed)
+                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
+                    FlatCombImg.addNoise(readout_noise)
 
                 FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
@@ -569,7 +595,7 @@ class Chip(FocalPlane):
                 datetime_obs = datetime.fromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
-                timestamp_obs += 5 * 60
+                timestamp_obs += 10 * 60
                 self.outputCal(
                     img=FlatCombImg,
                     ra_cen=ra_cen,
@@ -601,15 +627,15 @@ class Chip(FocalPlane):
             for i in range(NDark):
                 DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                     self.npix_x, self.npix_y, 
-                    overscan=overscan, bias_level=biaslevel, darkpsec=0.02, exptime=self.dark_exptime,
+                    overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time,
                     ncombine=1, read_noise=self.read_noise, 
                     gain=1, seed_bias=SeedBiasNonuni+self.chipID)
                 if config["cosmic_ray"].lower() == "y":
                     if config["cray_differ"].lower() == "y":
                         cr_map = effects.produceCR_Map(
                             xLen=self.npix_x, yLen=self.npix_y, 
-                            exTime=self.dark_exptime, 
-                            cr_pixelRatio=0.003*self.dark_exptime/self.exptime, 
+                            exTime=self.dark_exptime+0.5*self.readout_time, 
+                            cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), 
                             gain=self.gain, 
                             attachedSizes=self.attachedSizes,
                             seed=SeedCosmicRay+pointing_ID*30+self.chipID+2)
@@ -638,19 +664,35 @@ class Chip(FocalPlane):
                 # Non-Linearity for Dark
                 if config["non_linear"].lower() == "y":
                     print("  Applying Non-Linearity on the Dark image", flush=True)
-                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=1.e-7, beta2=1.e-10)
+                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0)
 
                 if config["cte_trail"].lower() == "y":
                     DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)
 
+                # Add Hot Pixels or/and Dead Pixels
+                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
+                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
+                DarkCombImg = effects.DefectivePixels(DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)
+
                 # Apply Bad lines 
                 if config["add_badcolumns"].lower() == "y":
                     DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID)
 
-                # Add Hot Pixels or/and Dead Pixels
-                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
-                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
-                DarkCombImg = effects.DefectivePixels(DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=self.bias_level)
+                # Add Bias level
+                if config["add_bias"].lower() == "y":
+                    print("  Adding Bias level and 16-channel non-uniformity")
+                    # img += float(config["bias_level"])
+                    DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, 
+                        bias_level=biaslevel, 
+                        nsecy = 2, nsecx=8, 
+                        seed=SeedBiasNonuni+self.chipID)
+
+                # Add Read-out Noise
+                if config["add_readout"].lower() == "y":
+                    seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID
+                    rng_readout = galsim.BaseDeviate(seed)
+                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
+                    DarkCombImg.addNoise(readout_noise)
 
                 DarkCombImg = effects.ApplyGainNonUniform16(
                     DarkCombImg, gain=self.gain, 
@@ -667,7 +709,7 @@ class Chip(FocalPlane):
                 datetime_obs = datetime.fromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
-                timestamp_obs += 5 * 60
+                timestamp_obs += 10 * 60
                 self.outputCal(
                     img=DarkCombImg,
                     ra_cen=ra_cen,
@@ -680,41 +722,6 @@ class Chip(FocalPlane):
                     output_dir=chip_output.subdir,
                     exptime=self.dark_exptime)
             del DarkCombImg
-
-        # Apply Nonlinearity on the chip image
-        if config["non_linear"].lower() == "y":
-            print("  Applying Non-Linearity on the chip image", flush=True)
-            img = effects.NonLinearity(GSImage=img, beta1=1.e-7, beta2=1.e-10)
-
-        # Apply CCD Saturation & Blooming
-        if config["saturbloom"].lower() == "y":
-            print("  Applying CCD Saturation & Blooming")
-            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)
-
-        # Apply CTE Effect
-        if config["cte_trail"].lower() == "y":
-            print("  Apply CTE Effect")
-            img = effects.CTE_Effect(GSImage=img, threshold=27)
-
-        # Apply Bad lines 
-        if config["add_badcolumns"].lower() == "y":
-            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID)
-
-        # Add Hot Pixels or/and Dead Pixels
-        rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
-        badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
-        img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=self.bias_level)
-        
-
-        # Apply Gain & Quantization
-        print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
-        img = effects.ApplyGainNonUniform16(
-            img, gain=self.gain, 
-            nsecy = 2, nsecx=8, 
-            seed=SeedGainNonuni+self.chipID)
-        img.array[img.array > 65535] = 65535
-        img.replaceNegative(replace_value=0)
-        img.quantize()
         # img = galsim.ImageUS(img)
 
         # # 16 output channel, with each a single image file

ObservationSim/Instrument/Chip/Effects.py

@@ -5,6 +5,7 @@ from numpy.core.fromnumeric import mean, size
 from numpy.random import Generator, PCG64
 import math
 from numba import jit
+from astropy import stats
 
 
 def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, widthb=8, read_noise=5):
@@ -31,7 +32,7 @@ def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8,
     return newimg
 
 
-def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=500):
+def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=0):
     # Also called bad pixels, including hot pixels and dead pixels
     # Hot Pixel > 20e-/s
     # Dead Pixel < 70%*Mean
@@ -68,10 +69,11 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=
     return GSImage
 
 
-def BadColumns(GSImage, seed=20210309, chipid=1):
+def BadColumns(GSImage, seed=20240309, chipid=1):
     # Set bad column values
     ysize,xsize = GSImage.array.shape
     subarr = GSImage.array[int(ysize*0.1):int(ysize*0.12), int(xsize*0.1):int(xsize*0.12)]
+    subarr = stats.sigma_clip(subarr, sigma=4, cenfunc='median', maxiters=3, masked=False)
     meanimg = np.median(subarr)
     stdimg = np.std(subarr)
     seed += chipid
@@ -83,12 +85,16 @@ def BadColumns(GSImage, seed=20210309, chipid=1):
     nbadsecA,nbadsecD = rgn.integers(low=1, high=5, size=2)
     collen = rgcollen.integers(low=int(ysize*0.1), high=int(ysize*0.7), size=(nbadsecA+nbadsecD)) 
     xposit = rgxpos.integers(low=int(xsize*0.05), high=int(xsize*0.95), size=(nbadsecA+nbadsecD))
-    signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1
-    dn = rgdn.integers(low=meanimg*0.4, high=meanimg*0.8, size=(nbadsecA+nbadsecD))*signs
+    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
     for badcoli in range(nbadsecA):
-        GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] += np.random.random((collen[badcoli],1))*stdimg*3+dn[badcoli]
+        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):
-        GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] += np.random.random((collen[badcoli+nbadsecA],1))*stdimg*3+dn[badcoli+nbadsecA]
+        GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli+nbadsecA],1)))+dn[badcoli+nbadsecA])
     return GSImage
 
 
@@ -223,7 +229,7 @@ def PRNU_Img(xsize, ysize, sigma=0.01, seed=202101):
     return prnuimg
 
 
-def NonLinearity(GSImage, beta1=1E-7, beta2=1E-10):
+def NonLinearity(GSImage, beta1=5E-7, beta2=0):
     NonLinear_f = lambda x, beta_1, beta_2: x - beta_1*x*x + beta_2*x*x*x
     GSImage.applyNonlinearity(NonLinear_f, beta1, beta2)
     return GSImage

ObservationSim/Instrument/Filter.py

@@ -73,5 +73,5 @@ class Filter(object):
     def getPhotonE(self):
         return photonEnergy(self.effective_wavelength)
 
-    def getSkyNoise(self, exptime, gain):
+    def getSkyNoise(self, exptime, gain=1.):
         return self.sky_background * exptime / gain

ObservationSim/ObservationSim.py

@@ -97,6 +97,7 @@ class Observation(object):
                 # Load SED
                 if obj.type == 'star':
                     normF = chip.normF_star
+                    #continue
                     try:
                         obj.load_SED(
                             survey_type=chip.survey_type, 
@@ -108,6 +109,7 @@ class Observation(object):
                         continue
                 elif obj.type == 'galaxy': # or obj.type == quasar
                     normF = chip.normF_galaxy
+                    # continue
                     obj.load_SED(
                         sed_path=sed_dir, 
                         survey_type=chip.survey_type, 
@@ -116,6 +118,7 @@ class Observation(object):
                         target_filt=filt)
                 elif obj.type == 'quasar':
                     normF = chip.normF_galaxy
+                    # continue
                     obj.load_SED(
                         sed_path=sed_dir, 
                         survey_type=chip.survey_type, 
@@ -200,7 +203,7 @@ class Observation(object):
 
         # Detector Effects
         # ===========================================================
-        chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map)
+        # chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map)
         # chip.img = chip.addEffects(config=self.config, img=chip.img, chip_output=chip_output, filt=filt, pointing_ID=pointing_ID)
 
         # whether to output zero, dark, flat calibration images.
@@ -214,7 +217,8 @@ class Observation(object):
             img_rot=img_rot,
             pointing_ID=pointing_ID,
             timestamp_obs=timestamp_obs,
-            pointing_type=pointing_type)
+            pointing_type=pointing_type,
+            sky_map=sky_map)
         
         if pointing_type == 'MS':
             datetime_obs = datetime.fromtimestamp(timestamp_obs)
@@ -309,13 +313,20 @@ class Observation(object):
                 sed_dir=self.path_dict["SED_dir"])
             print("finished running chip#%d..."%(chip.chipID), flush=True)
 
-    def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=[1621915200], pointing_type=['MS'], img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None):
+    def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None):
         comm = MPI.COMM_WORLD
         ind_thread = comm.Get_rank()
         num_thread = comm.Get_size()
 
         nchips_per_fp = len(self.chip_list)
 
+        # TEMP
+        if len(timestamp_obs) == 1:
+            timestamp_obs = np.tile(timestamp_obs, len(ra_cen))
+            pointing_type = np.tile(pointing_type, len(ra_cen))
+        
+        timestamp_obs = timestamp_obs[pRange]
+        pointing_type = pointing_type[pRange]
         ra_cen = ra_cen[pRange]
         dec_cen = dec_cen[pRange]
 
@@ -325,11 +336,6 @@ class Observation(object):
         else:
             pStart = 0
 
-        # TEMP
-        if len(timestamp_obs) == 1:
-            timestamp_obs *= len(pRange)
-            pointing_type *= len(pRange)
-
         for ipoint in range(len(ra_cen)):
             for ichip in range(nchips_per_fp):
                 i = ipoint*nchips_per_fp + ichip

ObservationSim/PSF/PSFInterp/PSFInterp.py

@@ -171,8 +171,8 @@ class PSFInterp(PSFModel):
             img = galsim.ImageF(imPSF, scale=pixSize)
             gsp = galsim.GSParams(folding_threshold=folding_threshold)
             self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
-            gaussPSF = galsim.Gaussian(sigma=0.04, gsparams=gsp)
-            self.psf    = galsim.Convolve(gaussPSF, self.psf)
+            # gaussPSF = galsim.Gaussian(sigma=0.04, gsparams=gsp)
+            # self.psf    = galsim.Convolve(gaussPSF, self.psf)
             return self.PSFspin(x=px/0.01, y=py/0.01)
         return imPSF
 

ObservationSim/Pointing.py

@@ -9,6 +9,8 @@ data_dir = "/data/simudata/CSSOSDataProductsSims/data/"
 config_file = os.path.join(work_dir, "ObservationSim.cfg")
 config = ReadConfig(config_file)
 
+exp_time = float(config['exp_time'])
+
 # Read Pointing list
 pointing_file = os.path.join(data_dir, "pointing10_20210202.dat")
 f = open(pointing_file, 'r')
@@ -37,12 +39,16 @@ pointing_type = ['CAL']*ncal + pointing_type
 t0 = datetime(2021, 5, 25, 12, 0, 0)
 t = datetime.timestamp(t0)
 timestamp_obs = []
+delta_t = 10 # Time elapsed between exposures (minutes)
 for i in range(len(pointing_type)):
     timestamp_obs.append(t)
     if pointing_type[i] == 'CAL':
-        t += 3 * 5 * 60 # 3 calibration exposure
+        t += 3 * delta_t * 60 # 3 calibration exposure
     elif pointing_type[i] == 'MS':
-        t += 5 * 60
+        t += delta_t * 60
+
+timestamp_obs = np.array(timestamp_obs)
+pointing_type = np.array(pointing_type)
 
 # Define the range of pointing list
 pRange = range(0, 2)
\ No newline at end of file

ObservationSim/run.pbs

@@ -7,7 +7,7 @@
 ##mpdboot -n 10 -f  ./mpd.hosts
 
 ##PBS -j oe
-#PBS -l nodes=comput110:ppn=80
+#PBS -l nodes=comput108:ppn=80
 #####PBS -q longq
 #PBS -q batch
 #PBS -u fangyuedong

ObservationSim/runExposure.py

@@ -20,4 +20,4 @@ gc.enable()
 # Testing run pointing list (using MPI)
 obs = Observation(work_dir=work_dir, data_dir=data_dir)
 # obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type)
-obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=150.)
+obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=exp_time)

test/ObservationSim.cfg

@@ -7,6 +7,7 @@
 # Observation setting
 date_obs      210525                  # Observation date [yymmdd]
 time_obs      120000                  # Observation time [hhmmss]
+exp_time      150.                    # Exposure time [seconds]
 ra_center     60.0                    # Telesscope pointing center [degree] (Default)
 dec_center   -40.0                     
 psf_rcont     0.15,0.8                # Radius of the 80% flux concentration (for a Gaussian PSF)
@@ -18,7 +19,7 @@ survey_type   Photometric             # Survey simulation option:
                                       # 'Photometric': only run sims for photometric chips
                                       # 'Spectroscopic': only run sims for spectroscopic chps
                                       # 'All': run sims for all chips
-psf_model     Interp                  # which PSF model to use:
+psf_model     Interp                  # Which PSF model to use:
                                       # 'Gauss': simple gaussin profile
                                       # 'Interp': Interpolated PSF from measured data 
 
@@ -26,8 +27,9 @@ psf_model     Interp                  # which PSF model to use:
 # Instrument effects
 #===========================================================
 field_dist       y                    # Y/N: Whether to add field distortion
-abs_back         y                    # Y: add sky background + sky noise & dark + dark noise 
-                                      # N: only add sky noise and dark noise
+add_back         y                    # Y/N: Whether to add sky background
+add_dark         y                    # Y/N: Whether to add dark noise
+add_readout      y                    # Y/N: Whether to add read-out (Gaussian) noise
 add_bias         y                    # Y/N: Whether to add bias-level to image
 shutter_effect   y                    # Y/N: Whether to apply shutter effect
 shutter_output   n                    # Y/N: Whether to export shutter effect 16-bit image (<=65535)
@@ -56,6 +58,7 @@ NDark            1                    # Number of frames to be exported
 NFlat            1                    # Number of frames to be exported
 dark_exptime     300                  # Exposure time for dark current frame
 flat_exptime     150                  # Exposure time for flat field frame
+readout_time     40                   # Time for read-out by each channel
 
 #===========================================================
 # Shear method
@@ -85,6 +88,7 @@ seed_mock            12345678         # Seed for generating random values of ra,
 seed_star            121212121        # Seed for generating random redshift, sed_type for stars
 seed_gal             212121212        # Seed for generating random redshift, sed_type for galaxies
 seed_Av              121212           # Seed for generating random intrinsic extinction
+seed_poisson         20210601         # Seed for Poisson noise
 seed_CR              20210317         # Seed for generating random cosmic ray map
 seed_flat            20210101         # Seed for generating random flat field
 seed_prnu            20210102         # Seed for photo-response non-uniformity
@@ -94,3 +98,4 @@ seed_biasNonUniform  20210203         # Seed for Bias nonuniformity
 seed_rnNonUniform    20210204         # Seed for ReadNoise nonuniformity
 seed_badcolumns      20240309         # Initial Seed for Bad Columns
 seed_defective       20210304         # Seed for Defective(bad) pixels
+seed_readout         20210601         # Seed for Read-out gaussian noise
\ No newline at end of file