led model

Catalog/C6_50sqdeg.py

@@ -343,7 +343,7 @@ class Catalog(CatalogBase):
                 input_time_str=time_str
             )
         for istars in range(nstars):
-            # # (TEST)
+            # (TEST)
             # if istars > 100:
             #     break
 

ObservationSim/Config/ChipOutput.py

 import os
 import logging
+import ObservationSim.Config._util as _util
 
 class ChipOutput(object):
-    def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='MS', pointing_ID='0', subdir="./", prefix=""):
+    def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='SCI', pointing_ID='0', subdir="./", prefix=""):
         self.focal_plane = focal_plane
         self.chip = chip
         self.filt = filt
@@ -20,12 +21,15 @@ class ChipOutput(object):
             self.dec_cen = config["obs_setting"]["dec_center"]
 
         self.chipLabel = focal_plane.getChipLabel(chip.chipID)
-        self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
+        obs_id = _util.get_obs_id(img_type=pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_num=pointing_ID)
+        # self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
+        self.cat_name = "MSC_%s_chip_%s_filt_%s"%(obs_id, self.chipLabel, filt.filter_type) + ".cat"
 
         self.subdir = subdir
 
         # Setup logger for each chip
-        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
+        # logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
+        logger_filename = "MSC_%s_chip_%s_filt_%s"%(obs_id, self.chipLabel, filt.filter_type) + ".log"
         self.logger = logging.getLogger()
         fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
         fh.setLevel(logging.DEBUG)
@@ -56,7 +60,7 @@ class ChipOutput(object):
         self.hdr += additional_column_names
     
     def create_output_file(self):
-        if self.pointing_type == 'MS':
+        if self.pointing_type == 'SCI':
             self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
             self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
             if not self.hdr.endswith("\n"):

ObservationSim/Config/Header/ImageHeader.py

@@ -17,6 +17,7 @@ from astropy.coordinates import SkyCoord
 from astropy.wcs.utils import fit_wcs_from_points
 from astropy.time import Time
 from astropy import wcs
+from ObservationSim.Config._util import get_obs_id, get_file_type
 
 from datetime import datetime
 # import socket
@@ -341,7 +342,7 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
 
 
 #TODO project_cycle is temporary， is not in header defined, delete in future
-def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, chip_name="01"):
+def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, run_counter=0, chip_name="01"):
 
     # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
 
@@ -376,18 +377,21 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     h_prim['OBJ_RA'] = ra
     h_prim['OBJ_DEC'] = dec
 
-    obs_type = {'SCI': '01', 'BIAS': '04', 'DARK': '08', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'DARKPF':'09'}
+    # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
 
-    OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
+    # # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
+    # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + str(int(run_counter)).rjust(2, '0') + pointNum.rjust(5,'0')
+    OBS_id = get_obs_id(img_type=im_type, project_cycle=project_cycle, run_counter=run_counter, pointing_num=pointNum)
 
     h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
     h_prim['OBSID'] = OBS_id
     # h_prim['TELFOCUS'] = 'f/14'
     h_prim['EXPTIME'] = exptime
     
-    # Define file types
-    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF','DARKPF':'DARK'}
-    h_prim['FILETYPE'] = file_type[im_type]
+    # # Define file types
+    # file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
+    # h_prim['FILETYPE'] = file_type[im_type]
+    h_prim['FILETYPE'] = get_file_type(img_type=im_type)
 
     co = coord.SkyCoord(ra, dec, unit='deg')
 
@@ -422,7 +426,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
     # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
     #     k - 1].rjust(2, '0') + '_L0_V01'
-    h_prim['FILENAME'] = 'CSST_MSC_MS_' + file_type[im_type] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
+    h_prim['FILENAME'] = 'CSST_MSC_MS_' + h_prim['FILETYPE'] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
 
 
     h_prim['POSI0_X'] = sat_pos[0]

ObservationSim/Config/Pointing.py

@@ -2,8 +2,10 @@ import numpy as np
 import galsim
 from astropy.time import Time
 
+import ObservationSim.Instrument._util as _util
+
 class Pointing(object):
-    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='MS'):
+    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='SCI'):
         self.id = id
         self.ra = ra
         self.dec = dec
@@ -14,9 +16,33 @@ class Pointing(object):
         self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
         self.exp_time = exp_time
         self.pointing_type = pointing_type
+        self.survey_field_type = 'WIDE'
         self.jdt = 0. 
     
-    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
+    def get_full_depth_exptime(self, filter_type):
+        if self.survey_field_type == 'WIDE':
+            if filter_type in _util.SPEC_FILTERS:
+                return 150. * 4
+            else:
+                if filter_type.lower() in ['nuv', 'y']:
+                    return 150. * 4
+                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
+                    return 150. * 2
+                else:
+                    return max(150., self.exp_time) # [TODO] for FGS
+        elif self.survey_field_type == 'DEEP':
+            if filter_type in _util.SPEC_FILTERS:
+                return 250. * 4 * 4
+            else:
+                if filter_type.lower() in ['nuv', 'y']:
+                    return 250. * 4 * 4
+                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
+                    return 250. * 2 * 4
+                else:
+                    return max(150., self.exp_time) # [TODO] for FGS
+
+
+    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCI'):
         self.id = id
         col_len = len(columns)
         self.ra = float(columns[0])
@@ -38,5 +64,9 @@ class Pointing(object):
             self.sat_vy = float(columns[16])
             self.sat_vz = float(columns[17])
             self.exp_time = float(columns[18])
+            is_deep = float(columns[19])
+            # [TODO] Can also define other survey types
+            if is_deep != -1.0:
+                self.survey_field_type = "DEEP"
         else:
             self.timestamp = t

ObservationSim/Instrument/Chip/Chip.py

@@ -3,6 +3,7 @@ import os
 import numpy as np
 import pickle
 import json
+import ObservationSim.Instrument._util as _util
 from astropy.table import Table
 from numpy.random import Generator, PCG64
 from astropy.io import fits
@@ -142,9 +143,9 @@ class Chip(FocalPlane):
         return rowID, colID
 
     def _getSurveyType(self):
-        if self.filter_type in ["GI", "GV", "GU"]:
+        if self.filter_type in _util.SPEC_FILTERS:
             return "spectroscopic"
-        elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
+        elif self.filter_type in _util.PHOT_FILTERS:
             return "photometric"
         # elif self.filter_type in ["FGS"]:
         #     return "FGS"
@@ -190,7 +191,7 @@ class Chip(FocalPlane):
     def getChipFilter(self, chipID=None):
         """Return the filter index and type for a given chip #(chipID)
         """
-        filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
+        filter_type_list = _util.ALL_FILTERS
         if chipID == None:
             chipID = self.chipID
 
@@ -228,7 +229,7 @@ class Chip(FocalPlane):
         for i in range(4):
             x = xcen + sign_x[i] * self.npix_x / 2.
             y = ycen + sign_y[i] * self.npix_y / 2.
-            x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
+            x, y = _util.rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
             xmin, xmax = min(xmin, x), max(xmax, x)
             ymin, ymax = min(ymin, y), max(ymax, y)
         return galsim.BoundsD(xmin, xmax, ymin, ymax)
@@ -267,7 +268,7 @@ class Chip(FocalPlane):
         noise = self.dark_noise * exptime + self.read_noise**2
         return noise
 
-    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, logger=None):
+    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, tel=None, logger=None):
         # Set random seeds
         SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
         SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
@@ -325,7 +326,7 @@ class Chip(FocalPlane):
         #     img.addNoise(poisson_noise)
 
         # Add cosmic-rays
-        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS':
+        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='SCI':
             chip_utils.log_info(msg="  Adding Cosmic-Ray", logger=self.logger)
             img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, 
                                                     seed=SeedCosmicRay+pointing_ID*30+self.chipID)
@@ -339,6 +340,8 @@ class Chip(FocalPlane):
                 pointing_ID=pointing_ID,
                 output_dir=chip_output.subdir,
                 exptime=exptime,
+                project_cycle=config["project_cycle"],
+                run_counter=config["run_counter"],
                 timestamp=timestamp_obs)
             del crmap_gsimg
 
@@ -358,7 +361,11 @@ class Chip(FocalPlane):
         #     img.addNoise(dark_noise)
 
         # Add dark current & Poisson noise
+        InputDark = False
         if config["ins_effects"]["add_dark"] == True:
+            if InputDark:
+                img = chip_utils.add_inputdark(img=img, chip=self, exptime=exptime)
+            else:
                 img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise)
         else:
             img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise, dark_noise=0.)
@@ -391,13 +398,23 @@ class Chip(FocalPlane):
             chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
             img = effects.CTE_Effect(GSImage=img, threshold=27)
 
+        ### prescan & overscan
+        if config["ins_effects"]["add_prescan"] == True:
+            img = chip_utils.AddPreScan(GSImage=img, pre1=27, pre2=4, over1=71, over2=80)
+        ### 1*16 output
+        if config["ins_effects"]["format_output"] == True:
+            img = chip_utils.formatOutput(GSImage=img)
+            self.nsecy = 1
+            self.nsecx = 16
+                    
+
         # Add Bias level
         if config["ins_effects"]["add_bias"] == True:
             chip_utils.log_info(msg="  Adding Bias level and 16-channel non-uniformity", logger=self.logger)
             if config["ins_effects"]["bias_16channel"] == True:
                 img = effects.AddBiasNonUniform16(img, 
                     bias_level=float(self.bias_level), 
-                    nsecy = 2, nsecx=8, 
+                    nsecy = self.nsecy, nsecx=self.nsecx, 
                     seed=SeedBiasNonuni+self.chipID,
                     logger=self.logger)
             elif config["ins_effects"]["bias_16channel"] == False:
@@ -415,7 +432,7 @@ class Chip(FocalPlane):
         if config["ins_effects"]["gain_16channel"] == True:
             img, self.gain_channel = effects.ApplyGainNonUniform16(
                 img, gain=self.gain, 
-                nsecy = 2, nsecx=8, 
+                nsecy = self.nsecy, nsecx=self.nsecx, 
                 seed=SeedGainNonuni+self.chipID,
                 logger=self.logger)
         elif config["ins_effects"]["gain_16channel"] == False:
@@ -490,6 +507,8 @@ class Chip(FocalPlane):
                     pointing_ID=pointing_ID,
                     output_dir=chip_output.subdir,
                     exptime=0.0,
+                    project_cycle=config["project_cycle"],
+                    run_counter=config["run_counter"],
                     timestamp=timestamp_obs)
             del BiasCombImg
 
@@ -590,6 +609,8 @@ class Chip(FocalPlane):
                     pointing_ID=pointing_ID,
                     output_dir=chip_output.subdir,
                     exptime=self.flat_exptime,
+                    project_cycle=config["project_cycle"],
+                    run_counter=config["run_counter"],
                     timestamp=timestamp_obs)
 
             del FlatCombImg, FlatSingle, prnu_img
@@ -643,6 +664,8 @@ class Chip(FocalPlane):
                         pointing_ID=pointing_ID,
                         output_dir=chip_output.subdir,
                         exptime=self.dark_exptime,
+                        project_cycle=config["project_cycle"],
+                        run_counter=config["run_counter"],
                         timestamp=timestamp_obs)
                     del crmap_gsimg
 
@@ -709,6 +732,8 @@ class Chip(FocalPlane):
                     pointing_ID=pointing_ID,
                     output_dir=chip_output.subdir,
                     exptime=self.dark_exptime,
+                    project_cycle=config["project_cycle"],
+                    run_counter=config["run_counter"],
                     timestamp = timestamp_obs)
             del DarkCombImg
         # img = galsim.ImageUS(img)

ObservationSim/Instrument/Chip/ChipUtils.py

@@ -53,7 +53,7 @@ def getChipSLSConf(chipID):
     if chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
     return confFile
 
-def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
+def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., project_cycle=6, run_counter=0, timestamp = 1621915200):
     datetime_obs = datetime.utcfromtimestamp(timestamp)
     date_obs = datetime_obs.strftime("%y%m%d")
     time_obs = datetime_obs.strftime("%H%M%S")
@@ -68,6 +68,8 @@ def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime
         time_obs=time_obs,
         im_type = im_type,
         exptime=exptime,
+        project_cycle=project_cycle,
+        run_counter=run_counter,
         chip_name=str(chip.chipID).rjust(2, '0')
         )
     h_ext = generateExtensionHeader(
@@ -91,7 +93,7 @@ def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime
         readoutTime = chip.readout_time)
     return h_prim, h_ext
 
-def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
+def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., project_cycle=6, run_counter=0, timestamp = 1621915200):
     h_prim, h_ext = generateHeader(
         chip=chip,
         ra_cen=ra_cen,
@@ -100,6 +102,8 @@ def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_
         im_type=im_type,
         pointing_ID=pointing_ID,
         exptime=exptime,
+        project_cycle=project_cycle,
+        run_counter=run_counter,
         timestamp=timestamp)
     hdu1 = fits.PrimaryHDU(header=h_prim)
     hdu1.add_checksum()
@@ -168,7 +172,8 @@ def get_poisson(seed=0, sky_level=0.):
     return rng_poisson, poisson_noise
 
 def get_base_img(img, read_noise, readout_time, dark_noise, exptime=150.):
-    base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time)
+    # base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time)
+    base_level = dark_noise*(exptime+0.5*readout_time)
     base_img = base_level * np.ones_like(img.array)
     return base_img
 
@@ -181,13 +186,14 @@ def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=Non
     base_img = get_base_img(img=img, read_noise=read_noise, readout_time=chip.readout_time, dark_noise=dark_noise, exptime=exptime)
     img += base_img
     img.addNoise(poisson_noise)
-    img -= read_noise**2
+    # img -= read_noise**2
     return img, base_img
 
 def add_brighter_fatter(img):
     #Inital dynamic lib
     try:
         with pkg_resources.files('ObservationSim.Instrument.Chip.lib_bf').joinpath("libmoduleBF.so") as lib_path:
+            print('--1', lib_path)
             lib_bf = ctypes.CDLL(lib_path)
     except AttributeError:
         with pkg_resources.path('ObservationSim.Instrument.Chip.lib_bf', "libmoduleBF.so") as lib_path:
@@ -210,3 +216,86 @@ def add_brighter_fatter(img):
     img.array[:, :] = np.reshape(arr_imc, [nx, ny])
     del arr_ima, arr_imc
     return img
+
+def add_inputdark(img, chip, exptime):
+    fname = "/share/home/weichengliang/CSST_git/test_new_sim/csst-simulation/ObservationSim/Instrument/data/dark/dark_1000s_example_0.fits"
+    hdu = fits.open(fname)
+    #ny, nx = img.array.shape
+    #inputdark = np.zeros([ny, nx])
+    img.array[:, :] += hdu[0].data/hdu[0].header['exptime']*exptime
+    hdu.close()
+    del inputdark
+    return img
+
+def AddPreScan(GSImage, pre1=27, pre2=4, over1=71, over2=80, nsecy = 2, nsecx=8):
+    img= GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+
+    imgt=np.zeros([int(nsecy*nsecx), int(ny/nsecy+pre2+over2), int(nx/nsecx+pre1+over1)])
+    for iy in range(nsecy):
+        for ix in range(nsecx):
+            if iy % 2 == 0:
+                tx = ix
+            else:
+                tx = (nsecx-1)-ix
+            ty = iy 
+            chunkidx = int(tx+ty*nsecx) #chunk1-[1,2,3,4], chunk2-[5,6,7,8], chunk3-[9,10,11,12], chunk4-[13,14,15,16]
+
+            imgtemp = np.zeros([int(ny/nsecy+pre2+over2), int(nx/nsecx+pre1+over1)])
+            if int(chunkidx/4) == 0:
+                imgtemp[pre2:-over2, pre1:-over1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp
+            if int(chunkidx/4) == 1:
+                imgtemp[pre2:-over2, over1:-pre1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[:, ::-1]
+            if int(chunkidx/4) == 2:
+                imgtemp[over2:-pre2, over1:-pre1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[::-1, ::-1]
+            if int(chunkidx/4) == 3:
+                imgtemp[over2:-pre2, pre1:-over1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[::-1, :]
+
+    imgtx1 = np.hstack(imgt[:nsecx:,       :, :])
+    imgtx2 = np.hstack(imgt[:(nsecx-1):-1, :, :])
+
+    newimg = galsim.Image(int(nx+(pre1+over1)*nsecx), int(ny+(pre2+over2)*nsecy), init_value=0)
+    newimg.array[:, :] = np.concatenate([imgtx1, imgtx2])
+
+    newimg.wcs = GSImage.wcs
+    return newimg
+
+def formatOutput(GSImage, nsecy = 2, nsecx=8):
+    img = GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+    
+    imgt = np.zeros([int(nsecx*nsecy), dy, dx])
+    for iy in range(nsecy):
+        for ix in range(nsecx):
+            if iy % 2 == 0:
+                tx = ix
+            else:
+                tx = (nsecx-1)-ix
+            ty = iy
+            chunkidx = int(tx+ty*nsecx)
+            if int(chunkidx/4) == 0:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 1:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 2:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 3:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+    
+    imgttx0 = np.hstack(imgt[ 0:4:,    :,    :])
+    imgttx1 = np.hstack(imgt[ 4:8:,    :, ::-1])
+    imgttx2 = np.hstack(imgt[8:12:, ::-1, ::-1])
+    imgttx3 = np.hstack(imgt[12:16:,::-1,    :])
+    
+    newimg = galsim.Image(int(dx*nsecx*nsecy), dy, init_value=0)
+    newimg.array[:, :] = np.hstack([imgttx0, imgttx1, imgttx2, imgttx3])
+    return newimg
+

ObservationSim/Instrument/Filter.py

@@ -4,7 +4,7 @@ import os
 import numpy as np
 import gc
 
-from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag
+import ObservationSim.Instrument._util as _util
 from ObservationSim.Instrument.FilterParam import FilterParam
 from ObservationSim.Straylight import Straylight
 
@@ -26,7 +26,7 @@ class Filter(object):
         self.survey_type = self._getSurveyType()
 
     def _getSurveyType(self):
-        if self.filter_type in ["GI", "GV", "GU"]:
+        if self.filter_type in _util.SPEC_FILTERS:
             return "spectroscopic"
         else:
             return "photometric"
@@ -106,20 +106,20 @@ class Filter(object):
         return bandpass_full, bandpass_sub_list
 
     def getPhotonE(self):
-        return photonEnergy(self.effective_wavelength)
+        return _util.photonEnergy(self.effective_wavelength)
 
     def getSkyNoise(self, exptime, gain=1.):
         return self.sky_background * exptime / gain
 
     def setFilterStrayLightPixel(self,jtime = 2460843., sat_pos = np.array([0,0,0]), pointing_radec = np.array([0,0]), sun_pos = np.array([0,0,0])):
         sl = Straylight(jtime=jtime, sat_pos=sat_pos, pointing_radec=pointing_radec,sun_pos=sun_pos)
-        if self.filter_type in ["GU","GV","GI"]:
+        if self.filter_type in _util.SPEC_FILTERS:
             s_pix, spec = sl.calculateStrayLightGrating(grating = self.filter_type.upper())
             if s_pix>0.8:
                 s_pix = 0.8
             self.sky_background = s_pix
             self.zodical_spec = spec
-        elif self.filter_type.lower() in ["nuv","u","g","r","i","z","y"]:
+        elif self.filter_type.lower() in [x.lower for x in _util.PHOT_FILTERS]:
             s_pix = sl.calculateStrayLightFilter(filter=self.filter_type.lower())
             if s_pix>1:
                 s_pix = 1
@@ -130,7 +130,7 @@ class Filter(object):
         gc.collect()
 
     def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None):
-        if self.filter_type in ["GI", "GV", "GU"]:
+        if self.filter_type in _util.SPEC_FILTERS:
             return
         if chip is not None:
             pix_scale = chip.pix_scale
@@ -144,6 +144,6 @@ class Filter(object):
             full_well = 90000
         throughput_file = self.filter_type.lower() + '_throughput.txt'
 
-        self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)
+        self.mag_limiting, self.mag_saturation = _util.calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)
 
         print("for filter %s: mag_limiting: %.3f, mag_saturation: %.3f"%(self.filter_type, self.mag_limiting, self.mag_saturation))

ObservationSim/Instrument/_util.py

@@ -14,6 +14,10 @@ VC_A = 2.99792458e+18  # speed of light: A/s
 VC_M = 2.99792458e+8   # speed of light: m/s
 H_PLANK = 6.626196e-27 # Plank constant: erg s
 
+ALL_FILTERS = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
+PHOT_FILTERS = ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']
+SPEC_FILTERS = ["GI", "GV", "GU"]
+
 def rotate_conterclockwise(x0, y0, x, y, angle):
 	"""
     Rotate a point counterclockwise by a given angle around a given origin.

ObservationSim/MockObject/FlatLED.py

@@ -98,7 +98,7 @@ class FlatLED(object):
             filt_bp = self.filt.filter_bandpass
             fil_eff = filt_bp.__call__(w_list / 10.)
             t_spec = np.trapz(f_spec*ccd_eff*fil_eff, w_list)
-            print(i, np.mean(unitFlatImg), t_spec, exp_t)
+            # print(i, np.mean(unitFlatImg), t_spec, exp_t)
             unitFlatImg = unitFlatImg * t_spec
 
             ledFlat = ledFlat+unitFlatImg*exp_t

ObservationSim/MockObject/MockObject.py

@@ -150,9 +150,10 @@ class MockObject(object):
             # Get PSF model
             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(nphotons)
-            star = galsim.Convolve(psf, star)
+            # 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)
             if np.sum(np.isnan(stamp.array)) > 0:

ObservationSim/ObservationSim.py

@@ -147,7 +147,7 @@ class Observation(object):
             sky_map = sky_map+filt.sky_background
             del flat_normal
 
-        if pointing.pointing_type == 'MS':
+        if pointing.pointing_type == 'SCI':
             # Load catalogues and templates
             self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt)
             chip_output.create_output_file()
@@ -156,7 +156,7 @@ class Observation(object):
             for ifilt in range(len(self.all_filter)):
                 temp_filter = self.all_filter[ifilt]
                 # Update the limiting magnitude using exposure time in pointing
-                temp_filter.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
+                temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip)
 
                 # Select cutting band filter for saturation/limiting magnitude
                 if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower():
@@ -197,14 +197,14 @@ class Observation(object):
             for j in range(self.nobj):
                 
                 # (DEBUG)
-                # if j >= 10:
-                #     break
+                if j >= 10:
+                    break
 
                 obj = self.cat.objs[j]
 
                 # (DEBUG)
-                if obj.getMagFilter(filt)>20:
-                    continue
+                # if obj.getMagFilter(filt)>20:
+                #     continue
 
                 # load and convert SED; also caculate object's magnitude in all CSST bands
                 try:
@@ -354,7 +354,7 @@ class Observation(object):
                 sky_map=sky_map, tel = self.tel,
                 logger=chip_output.logger)
             
-            if pointing.pointing_type == 'MS':
+            if pointing.pointing_type == 'SCI':
                 datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
@@ -371,6 +371,8 @@ class Observation(object):
                     im_type='SCI',
                     sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
                     sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
+                    project_cycle=self.config["project_cycle"],
+                    run_counter=self.config["run_counter"],
                     chip_name=str(chip.chipID).rjust(2, '0'))
                 h_ext = generateExtensionHeader(
                     chip=chip,
@@ -453,7 +455,6 @@ class Observation(object):
             sky_map=sky_map, tel=self.tel,
             logger=chip_output.logger)
 
-        if pointing.pointing_type == 'MS':
         datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
         date_obs = datetime_obs.strftime("%y%m%d")
         time_obs = datetime_obs.strftime("%H%M%S")

ObservationSim/_util.py

@@ -77,7 +77,7 @@ def generate_pointing_list(config, pointing_filename=None, data_dir=None):
             img_pa=config["obs_setting"]["image_rot"],
             timestamp=t,
             exp_time=exp_time,
-            pointing_type='MS'
+            pointing_type='SCI'
             )
         t += delta_t * 60.
         pointing_list.append(pointing)

config/config_50sqdeg.yaml

@@ -9,13 +9,13 @@
 # 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
-work_dir: "/share/home/fangyuedong/new_sim/workplace/"
+work_dir: "/share/home/weichengliang/CSST_git/test_new_sim/outputs/"
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "test_new_sim"
+run_name: "testRun0"
 
 # Whether to use MPI
 run_option:
-  use_mpi: NO
+  use_mpi: YES
   # 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
@@ -85,7 +85,7 @@ obs_setting:
   # if you just want to run default pointing:
   # - pointing_dir: null
   # - pointing_file: null
-  pointing_dir: "/share/home/fangyuedong/50sqdeg_pointings/"
+  pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/"
   pointing_file: "pointing_50_combined.dat"
 
   # Number of calibration pointings
@@ -163,6 +163,7 @@ ins_effects:
   add_dark:       ON  # Whether to add dark noise
   add_readout:    ON  # Whether to add read-out (Gaussian) noise
   add_bias:       ON  # Whether to add bias-level to images
+  add_prescan:    OFF
   bias_16channel: ON  # Whether to add different biases for 16 channels
   gain_16channel: ON  # Whether to make different gains for 16 channels
   shutter_effect: ON  # Whether to add shutter effect
@@ -177,6 +178,7 @@ ins_effects:
   add_hotpixels:  ON  # Whether to add hot pixels
   add_deadpixels: ON  # Whether to add dead(dark) pixels
   bright_fatter:  ON  # Whether to simulate Brighter-Fatter (also diffusion) effect
+  format_output:  OFF
 
   # Values: 
   # default values have been defined individually for each chip in:

config/config_C6.yaml

@@ -9,9 +9,12 @@
 # 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
-work_dir: "/share/home/fangyuedong/csst-simulation/workplace/"
+work_dir: "/share/home/fangyuedong/new_sim/workplace/"
+# work_dir: "/share/C6_new_sim_2sq"
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "profile_C6"
+run_name: "C6_new_sim_2sq_run1"
+project_cycle: 6
+run_counter: 1
 
 # Whether to use MPI
 run_option:
@@ -44,7 +47,7 @@ catalog_options:
     AGN_SED_WAVE: "wave_ross13.npy"
 
   # Only simulate stars?
-  star_only: NO
+  star_only: YES
 
   # Only simulate galaxies?
   galaxy_only: NO
@@ -112,7 +115,8 @@ obs_setting:
   cut_in_band: "z"
 
   # saturation magnitude margin
-  mag_sat_margin: -2.5
+  # mag_sat_margin: -2.5
+  mag_sat_margin: -15.
 
   # limiting magnitude margin
   mag_lim_margin: +1.0

config/config_C6_dev.yaml

@@ -11,7 +11,9 @@
 # OK to pass either way or both, as long as they are consistent
 work_dir: "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/"
 data_dir: "/Volumes/EAGET/C6_data/inputData/"
-run_name: "profile_C6"
+run_name: "C6_new_sim_2sq_run1"
+project_cycle: 6
+run_counter: 1
 
 # Whether to use MPI
 run_option:
@@ -44,7 +46,7 @@ catalog_options:
     AGN_SED_WAVE: "wave_ross13.npy"
 
   # Only simulate stars?
-  star_only: NO
+  star_only: YES
 
   # Only simulate galaxies?
   galaxy_only: NO
@@ -65,10 +67,10 @@ obs_setting:
   # "FGS": simulate FGS chips only (31-42)
   # "All": simulate full focal plane
   # "LED":
-  survey_type: "Photometric"
+  survey_type: "LED"
 
-  LED_TYPE: ["LED2","LED4","LED5"]
-  LED_TIME: [10.,10.,10.]
+  LED_TYPE: ["LED5"]
+  LED_TIME: [10.]
   
   # Exposure time [seconds]
   exp_time: 150.
@@ -98,7 +100,7 @@ obs_setting:
   # - 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: [80]
+  run_pointings: [0]
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
@@ -116,7 +118,8 @@ obs_setting:
   cut_in_band: "z"
 
   # saturation magnitude margin
-  mag_sat_margin: -2.5
+  # mag_sat_margin: -2.5
+  mag_sat_margin: -15.
 
   # limiting magnitude margin
   mag_lim_margin: +1.0
@@ -139,8 +142,7 @@ psf_setting:
   # path to PSF data
   # NOTE: only valid for "Interp" PSF
   psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
-  psf_sls_dir: "/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/"
-
+  psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/"
 ###############################################
 # Shear setting
 ###############################################
@@ -167,6 +169,7 @@ ins_effects:
   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
+  add_prescan:    OFF
   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
@@ -174,13 +177,14 @@ ins_effects:
   prnu_effect:    YES  # Whether to add PRNU effect
   non_linear:     YES  # Whether to add non-linearity
   cosmic_ray:     NO  # Whether to add cosmic-ray
-  cray_differ:    YES  # Whether to generate different cosmic ray maps CAL and MS output
+  cray_differ:    NO  # 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: NO  # Whether to add bad columns
-  add_hotpixels:  NO  # Whether to add hot pixels
-  add_deadpixels: NO  # Whether to add dead(dark) pixels
+  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
+  format_output:  OFF
 
   # Values: 
   # default values have been defined individually for each chip in:

profile_C6.sh

@@ -3,17 +3,22 @@
 date
 
 python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \
-    --config_file config_50sqdeg.yaml \
-    --catalog C6_50sqdeg \
+    --config_file config_C6.yaml \
+    --catalog C6_Catalog \
     -c /share/home/fangyuedong/new_sim/csst-simulation/config
+    
+    # --config_file config_test_new_sim.yaml \
+    # --catalog New_sim_Test \
+    # -c /share/home/fangyuedong/new_sim/csst-simulation/config
+    
+    # --config_file config_50sqdeg.yaml \
+    # --catalog C6_50sqdeg \
+    # -c /share/home/fangyuedong/new_sim/csst-simulation/config
+    
     # --config_file config_fgs.yaml \
     # --catalog FGS_Catalog \
     # -c /share/home/fangyuedong/csst-simulation/config
 
-    # --config_file config_C6.yaml \
-    # --catalog C6_Catalog \
-    # -c /share/home/fangyuedong/csst-simulation/config
-
     # --config_file test_fd_C6.yaml \
     # --catalog fd_test_C6 \
     # --config_file config_C6_test_wcs.yaml \

run_C6.pbs

@@ -8,16 +8,17 @@
 
 cd $PBS_O_WORKDIR
 
-NP=96
-hostfile=wcl-1,wcl-2
+NP=384
+hostfile=wcl-1,wcl-2,wcl-3,wcl-4,wcl-5,wcl-6
 date
 
 mpirun --oversubscribe -H $hostfile -np $NP python /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \
-    --config_file config_50sqdeg.yaml \
-    --catalog C6_50sqdeg \
+    --config_file config_C6.yaml \
+    --catalog C6_Catalog \
     -c /share/home/fangyuedong/new_sim/csst-simulation/config
-    # --config_file config_C6.yaml \
-    # --catalog C6_Catalog \
-    # -c /share/home/fangyuedong/csst-simulation/config
+
+    # --config_file config_50sqdeg.yaml \
+    # --catalog C6_50sqdeg \
+    # -c /share/home/fangyuedong/new_sim/csst-simulation/config
 
 

run_sim.py

@@ -58,6 +58,10 @@ def run_sim():
         config["obs_setting"]["mag_sat_margin"] = -2.5
     if "mag_lim_margin" not in config["obs_setting"]:
         config["obs_setting"]["mag_lim_margin"] = 1.0
+    if "project_cycle" not in config:
+        config["project_cycle"] = 6
+    if "run_counter" not in config:
+        config["run_counter"] = 0
 
     # Generate lists pointings based on the input pointing list (or default 
     # pointing RA, DEC) and "config["obs_setting"]["run_pointings"]".

setup.py

@@ -38,7 +38,7 @@ extensions = [
 
 df_module = [CTypes('ObservationSim.Instrument.Chip.lib_bf.libmoduleBF',
                     ['ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c', 'ObservationSim/Instrument/Chip/lib_bf/nrutil.c'],
-                    include_dirs=['ObservationSim/Instrument/Chip/lib_df_bf/', '/usr/include']
+                    include_dirs=['ObservationSim/Instrument/Chip/lib_bf/', '/usr/include']
 )]