issue17:fix led model memory bug; issue19:header time is not utc bug;...

issue17:fix led model memory bug; issue19:header time is not utc bug; issue18:header data-obs\darktime\exposure time fix by shaoli's chart;add new star catalog:C6_50sqdeg_ns.py

Catalog/C6_50sqdeg_ns.py

+import os
+import galsim
+import random
+import copy
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+from astropy.coordinates import spherical_to_cartesian
+from astropy.table import Table
+from scipy import interpolate
+from datetime import datetime
+
+from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+import astropy.constants as atcons
+
+import ctypes
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5']
+
+qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits']
+
+# star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5']
+star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)]
+
+star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5']
+
+class StarParm(ctypes.Structure):
+    _fields_ = [
+               ('logte',ctypes.c_float),
+               ('logg',ctypes.c_float),
+               ('Mass',ctypes.c_float),
+               ('Av', ctypes.c_float),
+               ('mu0', ctypes.c_float),
+               ('Z', ctypes.c_float)]
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+def get_agnsed_file(bundle_file_name):
+    return qsosed_file_list[bundle_file_list.index(bundle_file_name)]
+
+def get_star_cat(ra_pointing, dec_pointing):
+    pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg)
+    max_dist = 10
+    return_star_path = None
+    for star_file, center in zip(star_file_list, star_center_list):
+        center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg)
+        dist = pointing_c.separation(center_c).to(U.deg).value
+        if dist < max_dist:
+            return_star_path = star_file
+            max_dist = dist
+    return return_star_path
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = config["catalog_options"]["input_path"]["cat_dir"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            # Get the cloest star catalog file
+            star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec)
+            star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name)
+            self.star_path = os.path.join(self.cat_dir, star_path)
+            self.star_SED_path = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self._load_SED_lib_star()
+        
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = config["catalog_options"]["SED_templates_path"]["galaxy_SED"]
+            self._load_SED_lib_gals()
+            self.agn_seds = {}
+
+        if "AGN_SED" in config["catalog_options"]["SED_templates_path"] and not config["catalog_options"]["star_only"]:
+            self.AGN_SED_path = config["catalog_options"]["SED_templates_path"]["AGN_SED"]
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " av stellarmass dm teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_output_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    # def _load_SED_lib_star(self):
+    #     self.tempSED_star = h5.File(self.star_SED_path,'r')
+    def _load_SED_lib_star(self):
+        # self.tempSED_star = h5.File(self.star_SED_path,'r')
+        with pkg_resources.path('Catalog.data', 'starSpecInterp.so') as ddl_path:
+            self.starDDL = ctypes.CDLL(str(ddl_path))
+        self.starDDL.loadSpecLibs.argtypes=[ctypes.c_char_p, ctypes.c_char_p]
+        self.starDDL.loadExts.argtypes=[ctypes.c_char_p]
+        nwv = self.starDDL.loadSpecLibs(str.encode(os.path.join(self.star_SED_path,'file_BT-Settl_CSST_wl1000-24000_R1000.par')),str.encode(self.star_SED_path))
+        self.starDDL.loadExts(str.encode(os.path.join(self.star_SED_path,"Ext_odonnell94_R3.1_CSST_wl1000-24000_R1000.fits")))
+        self.star_spec_len = nwv
+
+    def _interp_star_sed(self, obj):
+        spec = (ctypes.c_float*self.star_spec_len)()
+        wave = (ctypes.c_float*self.star_spec_len)()
+        self.starDDL.interpSingleStar.argtypes=[ctypes.Structure, ctypes.POINTER(ctypes.c_float)]
+        # s=Star(obj.param['teff'], obj.param['grav''], obj.paramstar['mwmsc_mass'], obj.param['AV'], obj.param['DM'], obj.param['z_met'])
+        s=StarParm(obj.param['teff'], obj.param['logg'], obj.param['stellarMass'], obj.param['av'], obj.param['DM'], obj.param['feh'])
+
+        self.starDDL.interpSingleStar(s, spec, wave)
+
+        rv_c = obj.param['rv']/(atcons.c.value/1000.)
+        Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c))
+        wave_RV = wave*Doppler_factor
+        return wave_RV, spec
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0, agnsed_file=""):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=ngals,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        for igals in range(ngals):
+            # # (TEST)
+            # if igals > 100:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            
+            # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            if self.filt.filter_type == 'NUV':
+                param['mag_use_normal'] = gals['mag_csst_nuv'][igals]
+            else:
+                param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            param['g1'] = gals['shear'][igals][0]
+            param['g2'] = gals['shear'][igals][1]
+            param['kappa'] = gals['kappa'][igals]
+            param['e1'] = gals['ellipticity_true'][igals][0]
+            param['e2'] = gals['ellipticity_true'][igals][1]
+            
+            # For shape calculation
+            param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity(
+                                                                    e1=gals['ellipticity_true'][igals][0],
+                                                                    e2=gals['ellipticity_true'][igals][1],
+                                                                    rotation=self.rotation,
+                                                                    unit='radians')
+            # param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            # phi_e = cmath.phase(complex(param['e1'], param['e2']))
+            # param['e1'] = param['ell_total'] * np.cos(phi_e + 2*self.rotation)
+            # param['e2'] = param['ell_total'] * np.sin(phi_e + 2*self.rotation)
+            
+            param['e1_disk'] = param['e1']
+            param['e2_disk'] = param['e2']
+            param['e1_bulge'] = param['e1']
+            param['e2_bulge'] = param['e2']
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            param['size'] = gals['size'][igals]
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            # TEMP
+            self.ids += 1
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+
+            # Is this an Quasar?
+            param['qsoindex'] = gals['qsoindex'][igals]
+            if param['qsoindex'] == -1:
+                param['star'] = 0   # Galaxy
+                param['agnsed_file'] = ""
+                obj = Galaxy(param, logger=self.logger)
+            else:
+                param_qso = copy.deepcopy(param)
+                param_qso['star'] = 2   # Quasar
+                param_qso['agnsed_file'] = agnsed_file
+                # First add QSO model
+                obj = Quasar(param_qso, logger=self.logger)
+                # Need to deal with additional output columns
+                obj.additional_output_str = self.add_fmt%(0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+                                                        0, 0.)
+                self.objs.append(obj)
+                # Then add host galaxy model
+                param['star'] = 0   # Galaxy
+                param['agnsed_file'] = ""
+                obj = Galaxy(param, logger=self.logger)
+            
+            # Need to deal with additional output columns for (host) galaxy
+            obj.additional_output_str = self.add_fmt%(0., 0., 0., 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['RA'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["DEC"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=nstars,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["DEC"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_arr[istars]
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            param['mag_use_normal'] = stars['app_sdss_g'][istars]
+            self.ids += 1
+
+            param['id'] = '%06d'%(int(pix_id)) + '%08d'%(istars)
+            # param['sed_type'] = istars
+            # param['model_tag'] = ''
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['z_met'][istars]
+            param['stellarMass'] = stars['mass'][istars]
+
+
+            param['av'] = stars['AV'][istars]
+            param['DM'] = stars['DM'][istars]
+            # param['z_met'] = stars['z_met'][istars]
+
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            
+            try:
+                obj = Star(param, logger=self.logger)
+            except Exception as e:
+                print(e)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["av"], param['stellarMass'], param['DM'], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['star_catalog']
+            for pix in self.pix_list:
+                try:
+                    stars = star_cat[str(pix)]
+                    self._load_stars(stars, pix_id=pix)
+                    del stars
+                except Exception as e:
+                    self.logger.error(str(e))
+                    # print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    bundle_file = "galaxies_C6_bundle{:06}.h5".format(bundleID)
+                    file_path = os.path.join(self.galaxy_path, bundle_file)
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+
+                    # Get corresponding AGN SED file
+                    agnsed_file = get_agnsed_file(bundle_file)
+                    agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file)
+                    self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data
+
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file)
+
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            # _, wave, flux = tag_sed(
+            #     h5file=self.tempSED_star,
+            #     model_tag=obj.param['model_tag'],
+            #     teff=obj.param['teff'],
+            #     logg=obj.param['logg'],
+            #     feh=obj.param['feh']
+            # )
+            wave, flux = self._interp_star_sed(obj)
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.agn_seds[obj.agnsed_file][int(obj.qsoindex)] * 1e-17
+                flux[flux < 0] = 0.
+                wave = self.lamb_gal * (1.0 + obj.z)
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/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'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

ObservationSim/Config/Header/ImageHeader.py

@@ -19,7 +19,7 @@ from astropy.time import Time
 from astropy import wcs
 from ObservationSim.Config._util import get_obs_id, get_file_type
 
-from datetime import datetime
+from datetime import datetime, timezone
 # import socket
 import platform
 import toml
@@ -352,7 +352,10 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointing_id = '00000001', po
     # ccdnum = str(k)
 
     datetime_obs = datetime.utcfromtimestamp(time_pt)
+    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
+    # print(datetime_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5])
     datetime_obs = datetime.utcfromtimestamp(np.round(datetime_obs.timestamp(), 1))
+    # print(datetime_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5])
     # date_obs = datetime_obs.strftime("%y%m%d")
     # time_obs = datetime_obs.strftime("%H%M%S%f")[:-5]
 
@@ -551,30 +554,32 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
     h_ext['PIXSCAL1'] = pixel_scale
     h_ext['PIXSCAL2'] = pixel_scale
     h_ext['EXPTIME'] = exptime
-    h_ext['DARKTIME'] = exptime + readoutTime
+    h_ext['DARKTIME'] = exptime
 
     datetime_obs = datetime.utcfromtimestamp(timestamp)
+    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
     tstart = Time(datetime_obs)
     t_shutter_os = tstart
     t_shutter_oe = Time(tstart.mjd + t_shutter_open / 86400., format="mjd")
     t_shutter_co = Time(tstart.mjd + exptime / 86400., format="mjd")
     t_shutter_ce = Time(tstart.mjd + (exptime + t_shutter_close) / 86400., format="mjd")
-    t_shutter_os1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_os.unix).timestamp(), 1))
+    t_shutter_os1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_os.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['SHTOPEN0'] = t_shutter_os1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
-    t_shutter_oe1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_oe.unix).timestamp(), 1))
+    t_shutter_oe1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_oe.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['SHTOPEN1'] = t_shutter_oe1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
-    t_shutter_co1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_co.unix).timestamp(), 1))
+    t_shutter_co1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_co.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['SHTCLOS0'] = t_shutter_co1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
-    t_shutter_ce1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_ce.unix).timestamp(), 1))
+    t_shutter_ce1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_ce.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['SHTCLOS1'] = t_shutter_ce1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
     
+
     tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd")
     tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd")
     # tstart1=tstart.datetime.replace(microsecond=round(tstart.datetime.microsecond, -5))
-    tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart_read.unix).timestamp(), 1))
+    tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart_read.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['ROTIME0'] = tstart1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
     # tend_read1 = tend_read.datetime.replace(microsecond=round(tend_read.datetime.microsecond, -5))
-    tend_read1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tend_read.unix).timestamp(), 1))
+    tend_read1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tend_read.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
     h_ext['ROTIME1'] = tend_read1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]
     # h_ext['POS_ANG'] = pa
     header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size, 

ObservationSim/MockObject/CatalogBase.py

@@ -56,6 +56,8 @@ class CatalogBase(metaclass=ABCMeta):
             "teff":0.,
             "logg":0.,
             "feh":0.,
+            "DM":0.,
+            "stellarMass":1.,
             # C6 galaxies parameters
             "e1":0.,
             "e2":0.,

ObservationSim/MockObject/FlatLED.py

@@ -11,6 +11,7 @@ from scipy.interpolate import griddata
 from astropy.table import Table
 from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 from scipy import interpolate
+import gc
 
 from ObservationSim.MockObject.MockObject import MockObject
 # from ObservationSim.Straylight import calculateSkyMap_split_g
@@ -96,9 +97,48 @@ class FlatLED(object):
         U = griddata(X_, Z_, (
             M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)],
             N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]),
-                     method='cubic')
+                     method='linear')
         U = U/np.mean(U)
         flatImage = U*fluxLED[led_type]*1000
+        gc.collect()
+        return flatImage
+
+
+        ###
+    ### return LED flat, e/s
+    ###
+    def getLEDImage1(self, led_type='LED1'):
+        # cwave = cwaves[led_type]
+        flat = fits.open(os.path.join(self.flatDir, 'model_' + cwaves_name[led_type] + 'nm.fits'))
+        xlen = flat[0].header['NAXIS1']
+        ylen = 601
+
+        i = self.chip.rowID - 1
+        j = self.chip.colID - 1
+
+        x = np.linspace(0, self.chip.npix_x, int(xlen/6.))
+        y = np.linspace(0, self.chip.npix_y, int(ylen/5.))
+        xx, yy = np.meshgrid(x, y)
+
+        a1 = flat[0].data[int(ylen*i/5.):int(ylen*i/5.)+int(ylen/5.), int(xlen*j/6.):int(xlen*j/6.)+int(xlen/6.)]
+        # z = np.sin((xx+yy+xx**2+yy**2))
+        # fInterp = interp2d(xx, yy, z, kind='linear')
+
+        X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None]))
+        Z_ = a1.flatten()
+
+        n_x = np.arange(0, self.chip.npix_x , 1)
+        n_y = np.arange(0, self.chip.npix_y, 1)
+
+        M, N = np.meshgrid(n_x, n_y)
+
+        U = griddata(X_, Z_, (
+            M[0:self.chip.npix_y, 0:self.chip.npix_x],
+            N[0:self.chip.npix_y, 0:self.chip.npix_x ]),
+                     method='linear')
+        U = U/np.mean(U)
+        flatImage = U*fluxLED[led_type]*1000
+        gc.collect()
         return flatImage
 
     def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]):
@@ -114,7 +154,9 @@ class FlatLED(object):
         for i in np.arange(len(led_type_list)):
             led_type = led_type_list[i]
             exp_t = exp_t_list[i]
-            unitFlatImg = self.getLEDImage(led_type=led_type)
+            # unitFlatImg = self.getLEDImage(led_type=led_type)
+            unitFlatImg = self.getLEDImage1(led_type=led_type)
+            # print("---------------TEST mem:",np.mean(unitFlatImg))
             led_wave = cwaves[led_type]
             led_fwhm = cwaves_fwhm[led_type]
             led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm)
@@ -134,6 +176,7 @@ class FlatLED(object):
 
             ledStat = ledStat[0:int(led_type[3:])-1]+nledStat+ledStat[int(led_type[3:]):]
             ledTimes[int(led_type[3:])-1] = exp_t * 1000
+            gc.collect()
         return ledFlat, ledStat, ledTimes
 
     def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]):
@@ -150,7 +193,9 @@ class FlatLED(object):
         for i in np.arange(len(led_type_list)):
             led_type = led_type_list[i]
             exp_t = exp_t_list[i]
-            unitFlatImg = self.getLEDImage(led_type=led_type)
+            # unitFlatImg = self.getLEDImage(led_type=led_type)
+            unitFlatImg = self.getLEDImage1(led_type=led_type)
+            # print("---------------TEST mem:",np.mean(unitFlatImg))
             ledFlat_ = unitFlatImg*exp_t
             ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type]
             ledFlat_.astype(np.float32)

ObservationSim/sim_steps/add_LED_flat.py

@@ -2,6 +2,11 @@ import numpy as np
 from ObservationSim.MockObject import FlatLED
 import galsim
 
+from astropy.time import Time
+from datetime import datetime, timezone
+
+import gc
+
 def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
         
     if not hasattr(self, 'h_ext'):
@@ -18,6 +23,7 @@ def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
             self.updateHeaderInfo(header_flag='ext', keys = ['LEDSTAT'], values = [ledstat])
             self.updateHeaderInfo(header_flag='ext', keys = ['LEDT01','LEDT02','LEDT03','LEDT04','LEDT05','LEDT06','LEDT07','LEDT08','LEDT09','LEDT10','LEDT11','LEDT12','LEDT13','LEDT14'], values = letts)
         
+
     if obs_param["shutter_effect"] == True:
         pf_map = pf_map * chip.shutter_img
         pf_map = np.array(pf_map, dtype='float32')
@@ -26,4 +32,20 @@ def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
         self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','',''])
 
     chip.img = chip.img + pf_map
+
+    # renew header info
+    datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
+    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
+    t_obs = Time(datetime_obs)
+    
+    ##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光
+    t_obs_renew = Time(t_obs.mjd - (2.+0.5 +0.5) / 86400., format="mjd")
+
+    t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
+    self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
+
+    #dark time : 曝光时间+刷新后等带时间0.5s+点亮灯后0.5s+关闭快门时间1.5s+管快门后关灯前0.5+关灯后读出前等待0.5s
+    self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+0.5+1.5+0.5+0.5+pointing.exp_time])
+
+    gc.collect()
     return chip, filt, tel, pointing
\ No newline at end of file

ObservationSim/sim_steps/add_objects.py

@@ -7,6 +7,9 @@ import galsim
 from ObservationSim._util import get_shear_field
 from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS
 
+from astropy.time import Time
+from datetime import datetime, timezone
+
 def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
 
     # Get exposure time
@@ -59,6 +62,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
     # Get chip WCS
     if not hasattr(self, 'h_ext'):
         _, _ = self.prepare_headers(chip=chip, pointing=pointing)
+    
     chip_wcs = galsim.FitsWCS(header = self.h_ext)
     
     # Loop over objects
@@ -211,4 +215,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
     chip.img *= flat_normal
     del flat_normal
 
+
+    # renew header info
+    datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
+    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
+    t_obs = Time(datetime_obs)
+    
+    ##ccd刷新2s,等待0.5s，开始曝光
+    t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd")
+
+    t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
+    self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
+
+    #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
+    self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time])
+
     return chip, filt, tel, pointing
\ No newline at end of file

ObservationSim/sim_steps/add_sky_background.py

@@ -3,6 +3,9 @@ import galsim
 from ObservationSim.Straylight import calculateSkyMap_split_g
 from ObservationSim.Instrument import FilterParam
 
+from astropy.time import Time
+from datetime import datetime, timezone
+
 def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
     
     # Get exposure time
@@ -120,4 +123,20 @@ def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
             chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
         else:
             chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)
+
+    # renew header info
+    datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
+    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
+    t_obs = Time(datetime_obs)
+    
+    ##ccd刷新2s,等待0.5s，开始曝光
+    t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd")
+
+    t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
+    self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
+
+    #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
+    self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time])
+
+
     return chip, filt, tel, pointing
\ No newline at end of file

ObservationSim/sim_steps/readout_output.py

@@ -5,6 +5,9 @@ from astropy.io import fits
 from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
 from ObservationSim.Instrument.Chip import Effects
 
+from astropy.time import Time
+from datetime import datetime, timezone
+
 def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param):
     self.chip_output.Log_info("Apply pre/over-scan")
     chip.img = chip_utils.AddPreScan(GSImage=chip.img,
@@ -42,7 +45,17 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param)
 
     if not hasattr(self, 'h_ext'):
         _, _ = self.prepare_headers(chip=chip, pointing=pointing)
-        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','',''])
+        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1','EXPTIME'], values = [False,'','','','',0.0])
+        # renew header info
+        datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
+        datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
+        t_obs = Time(datetime_obs)
+    
+        ##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光
+        t_obs_renew = Time(t_obs.mjd - 2. / 86400., format="mjd")
+
+        t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
+        self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
     
     gains1 = list(chip.gain_channel[0:8])
     gains2 = list(chip.gain_channel[8:])

config/config_overall_newStar.yaml

+---
+###############################################
+#
+#  Configuration file for CSST simulation
+#           Overall settings
+#  CSST-Sim Group, 2024/01/08
+#
+###############################################
+
+# 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: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/"
+run_name: "QSO_50sqdeg_test"
+
+# Project cycle and run counter are used to name the outputs
+project_cycle: 9
+run_counter: 1
+
+# Run options
+run_option:
+  use_mpi: NO
+
+  # Output catalog only?
+  # If yes, no imaging simulation will be run. Only the catalogs
+  # of corresponding footprints will be generated.
+  out_cat_only: NO
+
+###############################################
+# Catalog setting
+###############################################
+# Configure the input catalog: options should be implemented
+# in the corresponding (user defined) 'Catalog' class
+catalog_options:
+  input_path:
+    cat_dir: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/"
+    star_cat: "star_catalog/"
+    # galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/"
+
+  SED_templates_path:
+    star_SED: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/star_catalog/"
+    # galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/"
+    # AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/"
+
+  # Only simulate stars?
+  star_only: YES
+
+  # Only simulate galaxies?
+  galaxy_only: NO
+
+###############################################
+# Observation setting
+###############################################
+obs_setting:
+  # (Optional) a file of point list 
+  # if you just want to run default pointing:
+  # - pointing_dir: null
+  # - pointing_file: null
+  pointing_file: "/nfsdata/share/CSSOSDataProductsSims/data/pointing_50_1_n.dat"
+
+  obs_config_file: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/csst-simulation/config/obs_config_SCI_WIDE_phot.yaml"
+
+  # Run specific pointing(s):
+  # - 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: [1]
+
+  # Whether to enable astrometric modeling
+  enable_astrometric_model: NO
+
+  # Cut by saturation magnitude in which band?
+  cut_in_band: "z"
+
+  # saturation magnitude margin
+  mag_sat_margin: -2.5
+  # mag_sat_margin: -15.
+
+  # limiting magnitude margin
+  mag_lim_margin: +1.0
+
+###############################################
+# PSF setting
+###############################################
+psf_setting:
+  
+  # Which PSF model to use:
+  # "Gauss": simple gaussian profile
+  # "Interp": Interpolated PSF from sampled ray-tracing data
+  psf_model: "Interp"
+
+  # PSF size [arcseconds]
+  # radius of 80% energy encircled
+  # NOTE: only valid for "Gauss" PSF
+  # psf_rcont: 0.15
+
+  # path to PSF data
+  # NOTE: only valid for "Interp" PSF
+  # PSF models for photometry survey simulation
+  psf_pho_dir: "/nfsdata/share/CSSOSDataProductsSims/data/psfCube1"
+  # PSF models for slitless spectrum survey simulation
+  psf_sls_dir: "/nfsdata/share/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/"
+
+###############################################
+# Shear setting
+###############################################
+
+shear_setting:
+  # Options to generate mock shear field:
+  # "constant": all galaxies are assigned a constant reduced shear
+  # "catalog": get shear values from catalog
+  shear_type: "constant"
+
+  # For constant shear field
+  reduced_g1: 0.
+  reduced_g2: 0.
+
+###############################################
+# Output options
+###############################################
+output_setting:
+  output_format:  "channels" # Whether to export as 16 channels (subimages) with pre- and over-scan ("image"/"channels")
+  shutter_output: NO # Whether to export shutter effect 16-bit image
+  prnu_output:    NO # Whether to export the PRNU (pixel-to-pixel flat-fielding) files
+
+###############################################
+# Random seeds
+###############################################
+random_seeds:
+  seed_poisson:         20210601  # Seed for Poisson noise
+  seed_CR:              20210317  # Seed for generating random cosmic ray maps
+  seed_flat:            20210101  # Seed for generating random flat fields
+  seed_prnu:            20210102  # Seed for photo-response non-uniformity
+  seed_gainNonUniform:  20210202  # Seed for gain nonuniformity
+  seed_biasNonUniform:  20210203  # Seed for bias nonuniformity
+  seed_rnNonUniform:    20210204  # Seed for readout-noise nonuniformity
+  seed_badcolumns:      20240309  # 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

setup.py

@@ -80,7 +80,7 @@ setup(name='CSSTSim',
         'ObservationSim.Instrument.data.throughputs': ['*.txt', '*.dat'], 
         'ObservationSim.Instrument.data.sls_conf': ['*.conf', '*.fits'],
         'ObservationSim.Instrument.data.flatCube': ['*.fits'],
-        'Catalog.data': ['*.fits'],
+        'Catalog.data': ['*.fits','*.so'],
         'ObservationSim.Config.Header':['*.header','*.lst'],
         'ObservationSim.Straylight.data': ['*.dat'],
         'ObservationSim.Straylight.data.sky': ['*.dat'],