add log for each treahds, add astrometry for pointings, add Catalog class for NGP fields

Catalog/C3Catalog.py

@@ -28,6 +28,11 @@ class C3Catalog(CatalogBase):
         self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
         self.seed_Av = config["random_seeds"]["seed_Av"]
 
+        if "logger" in kwargs:
+            self.logger = kwargs["logger"]
+        else:
+            self.logger = None
+
         with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
                 self.normF_star = Table.read(str(filter_path))
         with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path:
@@ -63,6 +68,10 @@ class C3Catalog(CatalogBase):
         dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
         vertices = spherical_to_cartesian(1., dec, ra)
         self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, 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):
@@ -169,10 +178,10 @@ class C3Catalog(CatalogBase):
             param['id'] = gals['galaxyID'][igals]
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation)
+                obj = Galaxy(param, self.rotation, logger=self.logger)
                 self.objs.append(obj)
             if param['star'] == 2:
-                obj = Quasar(param)
+                obj = Quasar(param, logger=self.logger)
                 self.objs.append(obj)
 
     def _load_stars(self, stars, pix_id=None):
@@ -230,7 +239,7 @@ class C3Catalog(CatalogBase):
             param['feh'] = stars['feh'][istars]
             param['z'] = 0.0
             param['star'] = 1   # Star
-            obj = Star(param)
+            obj = Star(param, logger=self.logger)
             self.objs.append(obj)
 
     def _load(self, **kwargs):
@@ -250,6 +259,9 @@ class C3Catalog(CatalogBase):
                 gals = gals_cat[str(pix)]
                 self._load_gals(gals, pix_id=pix)
                 del gals
+        if self.logger is not None:
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
             print("number of objects in catalog: ", len(self.objs))
         del self.avGal
 

Catalog/NGPCatalog.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+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 seds, sed_assign, extAv, tag_sed, getObservedSED
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+class NGPCatalog(CatalogBase):
+    def __init__(self, config, chip, pointing, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
+        self.seed_Av = config["random_seeds"]["seed_Av"]
+
+        if "logger" in kwargs:
+            self.logger = kwargs["logger"]
+        else:
+            self.logger = None
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+                self.normF_star = Table.read(str(filter_path))
+        with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path:
+                self.normF_galaxy = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"] and not config["run_option"]["galaxy_only"]:
+            star_file = config["input_path"]["star_cat"]
+            star_SED_file = config["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["input_path"] and config["input_path"]["galaxy_cat"] and not config["run_option"]["star_only"]:
+            galaxy_file = config["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+        if "rotateEll" in config["shear_setting"]:
+            self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
+        else:
+            self.rotation = 0.
+
+        self._get_healpix_list()
+        self._load()
+
+    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]))
+        vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, 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
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path)
+
+    def _load_gals(self, gals, pix_id=None):
+        ngals = len(gals['galaxyID'])
+        self.rng_sedGal = random.Random()
+        self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
+        self.ud = galsim.UniformDeviate(pix_id)
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra_true'][:]
+        dec_arr = gals['dec_true'][:]
+        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.fromtimestamp(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="J2015.5",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        for igals in range(ngals):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra_true'][igals]
+            param['dec_orig'] = gals['dec_true'][igals]
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
+            if param['mag_use_normal'] >= 26.5:
+                continue
+            param['z'] = gals['redshift_true'][igals]
+            param['model_tag'] = 'None'
+            param['gamma1'] = 0
+            param['gamma2'] = 0
+            param['kappa'] = 0
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+            # sersicB = gals['sersic_bulge'][igals]
+            hlrMajB = gals['size_bulge_true'][igals]
+            hlrMinB = gals['size_minor_bulge_true'][igals]
+            # sersicD = gals['sersic_disk'][igals]
+            hlrMajD = gals['size_disk_true'][igals]
+            hlrMinD = gals['size_minor_disk_true'][igals]
+            aGal = gals['size_true'][igals]
+            bGal = gals['size_minor_true'][igals]
+            param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
+            param['theta'] = gals['position_angle_true'][igals]
+            param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
+            param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
+            param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
+            param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
+            param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
+
+            # Assign each galaxy a template SED
+            param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
+            param['redden'] = self.tempRed_gal[param['sed_type']]
+            param['av'] = self.avGal[int(self.ud()*self.nav)]
+            if param['sed_type'] <= 5:
+                param['av'] = 0.0
+                param['redden'] = 0
+            param['star'] = 0   # Galaxy
+            if param['sed_type'] >= 29:
+                param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
+                param['star'] = 2 # Quasar
+
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = gals['galaxyID'][igals]
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+                self.objs.append(obj)
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+                self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # 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.fromtimestamp(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="J2015.5",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            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]
+            if param['mag_use_normal'] >= 26.5:
+                continue
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.nav = 15005
+        self.avGal = extAv(self.nav, seed=self.seed_Av)
+        self.objs = []
+        self.ids = 0
+        if "star_cat" in self.config["input_path"] and self.config["input_path"]["star_cat"] and not self.config["run_option"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['catalog']
+            for pix in self.pix_list:
+                stars = star_cat[str(pix)]
+                self._load_stars(stars, pix_id=pix)
+                del stars
+        if "galaxy_cat" in self.config["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
+            gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
+            for pix in self.pix_list:
+                gals = gals_cat[str(pix)]
+                self._load_gals(gals, pix_id=pix)
+                del gals
+        if self.logger is not None:
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+        del self.avGal
+
+
+    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']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            sed_data = getObservedSED(
+                sedCat=self.tempSed_gal[obj.sed_type],
+                redshift=obj.z,
+                av=obj.param["av"],
+                redden=obj.param["redden"]
+            )
+            wave, flux = sed_data[0], sed_data[1]
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        # lamb = np.arange(2500, 10001 + 0.5, 0.5)
+        lamb = np.arange(2400, 11001 + 0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photo/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        del wave
+        del flux
+        return sed

ObservationSim/Config/ChipOutput.py

 import os
+import logging
 
 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=""):
@@ -20,9 +21,22 @@ class ChipOutput(object):
         self.chipLabel = focal_plane.getChipLabel(chip.chipID)
         self.img_name =  prefix + exp_name%(self.chipLabel, filt.filter_type)
         
-        self.cat_name = 'MSC_' +  config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat"
+        # self.cat_name = 'MSC_' +  config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat"
+
+        self.cat_name = "MSC_%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat"
+
         self.subdir = subdir
 
+        # Setup logger for each chip
+        logger_filename = "MSC_%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), 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)
+        self.logger.setLevel(logging.DEBUG)
+        formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+        fh.setFormatter(formatter)
+        self.logger.addHandler(fh)
+
         hdr1  = "obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type "
         hdr2  = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 "
         hdr3  = "sed_type av redden "
@@ -36,10 +50,10 @@ class ChipOutput(object):
         self.hdr = hdr1 + hdr2 + hdr3 + hdr4
         self.fmt = fmt1 + fmt2 + fmt3 + fmt4
 
-        print("pointing_type = %s\n"%(pointing_type))
+        self.logger.info("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.logger.info("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

@@ -17,7 +17,7 @@ except ImportError:
     import importlib_resources as pkg_resources
 
 class Chip(FocalPlane):
-    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None):
+    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None):
         # Get focal plane (instance of paraent class) info
         # TODO: use chipID to config individual chip?
         super().__init__()
@@ -34,6 +34,8 @@ class Chip(FocalPlane):
         self.flat_exptime = float(config["ins_effects"]['flat_exptime'])
         self.readout_time = float(config["ins_effects"]['readout_time'])
 
+        self.logger = logger
+
         # A chip ID must be assigned
         self.chipID = int(chipID)
         self._getChipRowCol()
@@ -294,7 +296,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', sky_map=None, tel=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='MS', sky_map=None, tel=None, logger=None):
         SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
         SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
         SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"])
@@ -310,6 +312,7 @@ class Chip(FocalPlane):
             BoolDeadPix = True
         else:
             BoolDeadPix = False
+        self.logger = logger
 
         # Add sky background
         if sky_map is None:
@@ -324,6 +327,11 @@ class Chip(FocalPlane):
 
         # Apply flat-field large scale structure for one chip
         if config["ins_effects"]["flat_fielding"] == True:
+            if self.logger is not None:
+                self.logger.info("  Creating and applying Flat-Fielding")
+                msg = str(img.bounds)
+                self.logger.info(msg)
+            else:
                 print("  Creating and applying Flat-Fielding", flush=True)
                 print(img.bounds, flush=True)
             flat_img = effects.MakeFlatSmooth(
@@ -338,6 +346,9 @@ class Chip(FocalPlane):
 
         # Apply Shutter-effect for one chip
         if config["ins_effects"]["shutter_effect"] == True:
+            if self.logger is not None:
+                self.logger.info("  Apply shutter effect")
+            else:
                 print("  Apply shutter effect", flush=True)
             shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3)    # shutter effect normalized image for this chip
             if self.survey_type == "photometric":
@@ -356,6 +367,9 @@ class Chip(FocalPlane):
 
         # Add cosmic-rays
         if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS':
+            if self.logger is not None:
+                self.logger.info(("  Adding Cosmic-Ray"))
+            else:
                 print("  Adding Cosmic-Ray", flush=True)
             cr_map, cr_event_num = effects.produceCR_Map(
                 xLen=self.npix_x, yLen=self.npix_y, 
@@ -389,6 +403,9 @@ class Chip(FocalPlane):
 
         # Apply PRNU effect and output PRNU flat file:
         if config["ins_effects"]["prnu_effect"] == True:
+            if self.logger is not None:
+                self.logger.info("  Applying PRNU effect")
+            else:
                 print("  Applying PRNU effect", flush=True)
             prnu_img = effects.PRNU_Img(
                 xsize=self.npix_x, 
@@ -413,28 +430,41 @@ class Chip(FocalPlane):
 
         # Apply Bad lines 
         if config["ins_effects"]["add_badcolumns"] == True:
-            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID)
+            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
 
         # Add Bias level
         if config["ins_effects"]["add_bias"] == True:
+            if self.logger is not None:
+                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
+            else:
                 print("  Adding Bias level and 16-channel non-uniformity")
             img = effects.AddBiasNonUniform16(img, 
                 bias_level=float(config["ins_effects"]["bias_level"]), 
                 nsecy = 2, nsecx=8, 
-                seed=SeedBiasNonuni+self.chipID)
+                seed=SeedBiasNonuni+self.chipID,
+                logger=self.logger)
 
         # Apply Nonlinearity on the chip image
         if config["ins_effects"]["non_linear"] == True:
+            if self.logger is not None:
+                self.logger.info("  Applying Non-Linearity on the chip image")
+            else:
                 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["ins_effects"]["saturbloom"] == True:
+            if self.logger is not None:
+                self.logger.info("  Applying CCD Saturation & Blooming")
+            else:
                 print("  Applying CCD Saturation & Blooming")
             img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)
 
         # Apply CTE Effect
         if config["ins_effects"]["cte_trail"] == True:
+            if self.logger is not None:
+                self.logger.info("  Apply CTE Effect")
+            else:
                 print("  Apply CTE Effect")
             img = effects.CTE_Effect(GSImage=img, threshold=27)
 
@@ -447,11 +477,15 @@ class Chip(FocalPlane):
 
 
         # Apply Gain & Quantization
+        if self.logger is not None:
+            self.logger.info("  Applying Gain (and 16 channel non-uniformity) & Quantization")
+        else:
             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)
+            seed=SeedGainNonuni+self.chipID,
+            logger=self.logger)
         img.array[img.array > 65535] = 65535
         img.replaceNegative(replace_value=0)
         img.quantize()
@@ -461,6 +495,9 @@ class Chip(FocalPlane):
         ######################################################################################
         # Bias output
         if config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
+            if self.logger is not None:
+                self.logger.info("  Output N frame Bias files")
+            else:
                 print("  Output N frame Bias files", flush=True)
             NBias = int(config["ins_effects"]["NBias"])
             for i in range(NBias):
@@ -468,7 +505,8 @@ class Chip(FocalPlane):
                     self.npix_x, self.npix_y, 
                     bias_level=float(config["ins_effects"]["bias_level"]), 
                     ncombine=1, read_noise=self.read_noise, gain=1,
-                    seed=SeedBiasNonuni+self.chipID)
+                    seed=SeedBiasNonuni+self.chipID,
+                    logger=self.logger)
                 if config["ins_effects"]["cosmic_ray"] == True:
                     if config["ins_effects"]["cray_differ"] == True:
                         cr_map, cr_event_num = effects.produceCR_Map(
@@ -484,16 +522,20 @@ class Chip(FocalPlane):
 
                 # Non-Linearity for Bias
                 if config["ins_effects"]["non_linear"] == True:
+                    if self.logger is not None:
+                        self.logger.info("  Applying Non-Linearity on the Bias image")
+                    else:
                         print("  Applying Non-Linearity on the Bias image", flush=True)
                     BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0)
 
                 # Apply Bad lines 
                 if config["ins_effects"]["add_badcolumns"] == True:
-                    BiasCombImg = effects.BadColumns(BiasCombImg-float(config["ins_effects"]["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID) + float(config["ins_effects"]["bias_level"])-5
+                    BiasCombImg = effects.BadColumns(BiasCombImg-float(config["ins_effects"]["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(config["ins_effects"]["bias_level"])-5
 
                 BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
-                    seed=SeedGainNonuni+self.chipID)
+                    seed=SeedGainNonuni+self.chipID,
+                    logger=self.logger)
                 # BiasCombImg = effects.AddOverscan(
                 #     BiasCombImg, 
                 #     overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain, 
@@ -521,6 +563,9 @@ class Chip(FocalPlane):
 
         # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
         if config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
+            if self.logger is not None:
+                self.logger.info("  Output N frame Flat-Field files")
+            else:
                 print("  Output N frame Flat-Field files", flush=True)
             NFlat = int(config["ins_effects"]["NFlat"])
             if config["ins_effects"]["add_bias"] == True:
@@ -539,7 +584,8 @@ class Chip(FocalPlane):
                     gain=1, 
                     overscan=overscan, 
                     biaslevel=0,
-                    seed_bias=SeedDefective+self.chipID
+                    seed_bias=SeedDefective+self.chipID,
+                    logger=self.logger
                     )
                 if config["ins_effects"]["cosmic_ray"] == True:
                     if config["ins_effects"]["cray_differ"] == True:
@@ -555,6 +601,9 @@ class Chip(FocalPlane):
                     del cr_map
 
                 if config["ins_effects"]["non_linear"] == True:
+                    if self.logger is not None:
+                        self.logger.info("  Applying Non-Linearity on the Flat image")
+                    else:
                         print("  Applying Non-Linearity on the Flat image", flush=True)
                     FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0)
 
@@ -568,16 +617,20 @@ class Chip(FocalPlane):
 
                 # Apply Bad lines 
                 if config["ins_effects"]["add_badcolumns"] == True:
-                    FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID)
+                    FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
 
                 # Add Bias level
                 if config["ins_effects"]["add_bias"] == True:
+                    if self.logger is not None:
+                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
+                    else:
                         print("  Adding Bias level and 16-channel non-uniformity")
                     # img += float(config["ins_effects"]["bias_level"])
                     FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, 
                         bias_level=biaslevel, 
                         nsecy = 2, nsecx=8, 
-                        seed=SeedBiasNonuni+self.chipID)
+                        seed=SeedBiasNonuni+self.chipID,
+                        logger=self.logger)
                 
                 # Add Read-out Noise
                 if config["ins_effects"]["add_readout"] == True:
@@ -588,7 +641,8 @@ class Chip(FocalPlane):
 
                 FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
-                    seed=SeedGainNonuni+self.chipID)
+                    seed=SeedGainNonuni+self.chipID,
+                    logger=self.logger)
                 # FlatCombImg = effects.AddOverscan(FlatCombImg, overscan=overscan, gain=self.gain, widthl=27, widthr=27, widtht=8, widthb=8)
                 FlatCombImg.replaceNegative(replace_value=0)
                 FlatCombImg.quantize()
@@ -618,6 +672,9 @@ class Chip(FocalPlane):
 
         # Export Dark current images
         if config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
+            if self.logger is not None:
+                self.logger.info("  Output N frame Dark Current files")
+            else:
                 print("  Output N frame Dark Current files", flush=True)
             NDark = int(config["ins_effects"]["NDark"])
             if config["ins_effects"]["add_bias"] == True:
@@ -631,7 +688,8 @@ class Chip(FocalPlane):
                     self.npix_x, self.npix_y, 
                     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)
+                    gain=1, seed_bias=SeedBiasNonuni+self.chipID,
+                    logger=self.logger)
                 if config["ins_effects"]["cosmic_ray"] == True:
                     if config["ins_effects"]["cray_differ"] == True:
                         cr_map, cr_event_num = effects.produceCR_Map(
@@ -665,6 +723,9 @@ class Chip(FocalPlane):
 
                 # Non-Linearity for Dark
                 if config["ins_effects"]["non_linear"] == True:
+                    if self.logger is not None:
+                        self.logger.info("  Applying Non-Linearity on the Dark image")
+                    else:
                         print("  Applying Non-Linearity on the Dark image", flush=True)
                     DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0)
 
@@ -678,16 +739,20 @@ class Chip(FocalPlane):
 
                 # Apply Bad lines 
                 if config["ins_effects"]["add_badcolumns"] == True:
-                    DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID)
+                    DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
 
                 # Add Bias level
                 if config["ins_effects"]["add_bias"] == True:
+                    if self.logger is not None:
+                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
+                    else:
                         print("  Adding Bias level and 16-channel non-uniformity")
                     # img += float(config["ins_effects"]["bias_level"])
                     DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, 
                         bias_level=biaslevel, 
                         nsecy = 2, nsecx=8, 
-                        seed=SeedBiasNonuni+self.chipID)
+                        seed=SeedBiasNonuni+self.chipID,
+                        logger=self.logger)
 
                 # Add Read-out Noise
                 if config["ins_effects"]["add_readout"] == True:
@@ -699,7 +764,8 @@ class Chip(FocalPlane):
                 DarkCombImg = effects.ApplyGainNonUniform16(
                     DarkCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
-                    seed=SeedGainNonuni+self.chipID)
+                    seed=SeedGainNonuni+self.chipID,
+                    logger=self.logger)
                 # DarkCombImg = effects.AddOverscan(
                 #     DarkCombImg, 
                 #     overscan=overscan, gain=self.gain, 

ObservationSim/Instrument/Chip/Effects.py

@@ -69,7 +69,7 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=
     return GSImage
 
 
-def BadColumns(GSImage, seed=20240309, chipid=1):
+def BadColumns(GSImage, seed=20240309, chipid=1, logger=None):
     # 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)]
@@ -85,6 +85,9 @@ def BadColumns(GSImage, seed=20240309, 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))
+    if logger is not None:
+        logger.info(xposit+1)
+    else:
         print(xposit+1)
     # signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1
     # if meanimg>0:
@@ -98,7 +101,7 @@ def BadColumns(GSImage, seed=20240309, chipid=1):
     return GSImage
 
 
-def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=202102):
+def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=202102, logger=None):
     # Generate Bias and its non-uniformity, and add the 16 bias values to the GS-Image
     rg = Generator(PCG64(int(seed)))
     Random16 = (rg.random(nsecy*nsecx)-0.5)*20
@@ -106,6 +109,10 @@ def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=2021
         BiasLevel = np.zeros((nsecy,nsecx))
     elif bias_level>0:
         BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level
+    if logger is not None:
+        msg = str(" Biases of 16 channels: " + str(BiasLevel))
+        logger.info(msg)
+    else:
         print(" Biases of 16 channels:\n",BiasLevel)
     arrshape = GSImage.array.shape
     secsize_x = int(arrshape[1]/nsecx)
@@ -116,14 +123,15 @@ def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=2021
     return GSImage
 
 
-def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain=1, seed=202102):
+def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain=1, seed=202102, logger=None):
     # Start with 0 value bias GS-Image
     ncombine=int(ncombine)
     BiasSngImg0 = galsim.Image(npix_x, npix_y, init_value=0)
     BiasSngImg = AddBiasNonUniform16(BiasSngImg0, 
                 bias_level=bias_level, 
                 nsecy = 2, nsecx=8, 
-                seed=int(seed))
+                seed=int(seed),
+                logger=logger)
     BiasCombImg = BiasSngImg*ncombine
     rng = galsim.UniformDeviate()
     NoiseBias = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5)
@@ -139,11 +147,15 @@ def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain
     return BiasCombImg, BiasTag
 
 
-def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
+def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None):
     # Generate Gain non-uniformity, and multipy the different factors (mean~1 with sigma~1%) to the GS-Image
     rg = Generator(PCG64(int(seed)))
     Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1   # sigma~1%
     Gain16 = Random16.reshape((nsecy,nsecx))/gain
+    if logger is not None:
+        msg = str("Gain of 16 channels: " + str(Gain16))
+        logger.info(msg)
+    else:
         print("Gain of 16 channels: ",Gain16)
     arrshape = GSImage.array.shape
     secsize_x = int(arrshape[1]/nsecx)
@@ -154,11 +166,15 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
     return GSImage
 
 
-def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
+def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None):
     # Generate Gain non-uniformity, and multipy the different factors (mean~1 with sigma~1%) to the GS-Image
     rg = Generator(PCG64(int(seed)))
     Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1   # sigma~1%
     Gain16 = Random16.reshape((nsecy,nsecx))/gain
+    if logger is not None:
+        msg = str(seed-20210202, "Gains of 16 channels: " + str(Gain16))
+        logger.info(msg)
+    else:
         print(seed-20210202, "Gains of 16 channels:\n", Gain16)
     # arrshape = GSImage.array.shape
     # secsize_x = int(arrshape[1]/nsecx)
@@ -186,7 +202,7 @@ def MakeFlatSmooth(GSBounds, seed):
     return FlatImg
 
 
-def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=500, biaslevel=500, seed_bias=20210311):
+def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=500, biaslevel=500, seed_bias=20210311, logger=None):
     ncombine=int(ncombine)
     FlatCombImg = flat_single_image*ncombine
     rng = galsim.UniformDeviate()
@@ -200,7 +216,8 @@ def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=
             FlatCombImg, 
             bias_level=biaslevel, 
             nsecy=2, nsecx=8, 
-            seed=seed_bias)
+            seed=seed_bias,
+            logger=logger)
     if ncombine == 1:
         FlatTag = 'Single'
         pass
@@ -212,7 +229,7 @@ def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=
     return FlatCombImg, FlatTag
 
 
-def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102, darkpsec=0.02, exptime=150, ncombine=10, read_noise=5, gain=1):
+def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102, darkpsec=0.02, exptime=150, ncombine=10, read_noise=5, gain=1, logger=None):
     ncombine=int(ncombine)
     darkpix = darkpsec*exptime
     DarkSngImg = galsim.Image(npix_x, npix_y, init_value=darkpix)
@@ -227,7 +244,8 @@ def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102
             DarkCombImg, 
             bias_level=bias_level, 
             nsecy = 2, nsecx=8, 
-            seed=int(seed_bias))
+            seed=int(seed_bias),
+            logger=logger)
     if ncombine == 1:
         DarkTag = 'Single'
         pass

ObservationSim/MockObject/Galaxy.py

@@ -10,8 +10,8 @@ from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 from ObservationSim.MockObject.MockObject import MockObject
 
 class Galaxy(MockObject):
-    def __init__(self, param, rotation=None):
-        super().__init__(param)
+    def __init__(self, param, rotation=None, logger=None):
+        super().__init__(param, logger=logger)
         self.thetaR = self.param["theta"]
         self.bfrac = self.param["bfrac"]
         self.hlr_disk = self.param["hlr_disk"]

ObservationSim/MockObject/MockObject.py

@@ -8,7 +8,7 @@ from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFacto
 from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 
 class MockObject(object):
-    def __init__(self, param):
+    def __init__(self, param, logger=None):
         self.param = param
 
         if self.param["star"] == 0:
@@ -45,6 +45,8 @@ class MockObject(object):
         self.e1_disk, self.e2_disk = 0., 0.
         self.e1_bulge, self.e2_bulge = 0., 0.
 
+        self.logger = logger
+
     def getMagFilter(self, filt):
         if filt.filter_type in ["GI", "GV", "GU"]:
             return self.param["mag_use_normal"]

ObservationSim/MockObject/Quasar.py

@@ -9,8 +9,8 @@ from ObservationSim.MockObject.MockObject import MockObject
 from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
 
 class Quasar(MockObject):
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, logger=None):
+        super().__init__(param, logger=logger)
 
     def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
         if survey_type == "photometric":

ObservationSim/MockObject/Star.py

@@ -9,8 +9,8 @@ from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFacto
 from ObservationSim.MockObject.MockObject import MockObject
 
 class Star(MockObject):
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, logger=None):
+        super().__init__(param, logger=logger)
 
     def unload_SED(self):
         """(Test) free up SED memory

ObservationSim/ObservationSim.py

@@ -14,6 +14,7 @@ from ObservationSim.Instrument.Chip import Effects
 from ObservationSim.MockObject import calculateSkyMap_split_g
 from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
 from ObservationSim._util import get_shear_field, makeSubDir_PointingList
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
 
 class Observation(object):
     def __init__(self, config, Catalog, work_dir=None, data_dir=None):
@@ -56,30 +57,68 @@ class Observation(object):
 
     def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, shear_cat_file=None, cat_dir=None, sed_dir=None):
 
-        print(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
-        print("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
-        print("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat())
-        print("Exposure time: %f" % pointing.exp_time)
-        print("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
-        print("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
-        print("Position Angle: %f" % pointing.img_pa.deg)
-        print('Chip : %d' % chip.chipID)
-        print(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
+        # print(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
+        # print("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
+        # print("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat())
+        # print("Exposure time: %f" % pointing.exp_time)
+        # print("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
+        # print("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
+        # print("Position Angle: %f" % pointing.img_pa.deg)
+        # print('Chip : %d' % chip.chipID)
+        # print(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
+        chip_output.logger.info(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
+        chip_output.logger.info("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
+        chip_output.logger.info("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat())
+        chip_output.logger.info("Exposure time: %f" % pointing.exp_time)
+        chip_output.logger.info("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
+        chip_output.logger.info("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
+        chip_output.logger.info("Position Angle: %f" % pointing.img_pa.deg)
+        chip_output.logger.info('Chip : %d' % chip.chipID)
+        chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
 
         if self.config["psf_setting"]["psf_model"] == "Gauss":
             psf_model = PSFGauss(chip=chip)
         elif self.config["psf_setting"]["psf_model"] == "Interp":
             psf_model = PSFInterp(chip=chip, PSF_data_file=self.path_dict["psf_dir"])
         else:
-            print("unrecognized PSF model type!!", flush=True)
+            # print("unrecognized PSF model type!!", flush=True)
+            chip_output.logger.error("unrecognized PSF model type!!", flush=True)
 
         # Get (extra) shear fields
         if shear_cat_file is not None:
             self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config, shear_cat_file=shear_cat_file)
 
+        # Apply astrometric simulation for pointing
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            dt = datetime.fromtimestamp(pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_cen, dec_cen = on_orbit_obs_position(
+                input_ra_list=[pointing.ra],
+                input_dec_list=[pointing.dec],
+                input_pmra_list=[0.],
+                input_pmdec_list=[0.],
+                input_rv_list=[0.],
+                input_parallax_list=[1e-9],
+                input_nstars=1,
+                input_x=pointing.sat_x,
+                input_y=pointing.sat_y,
+                input_z=pointing.sat_z,
+                input_vx=pointing.sat_vx,
+                input_vy=pointing.sat_vy,
+                input_vz=pointing.sat_vz,
+                input_epoch="J2015.5",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+            ra_cen, dec_cen = ra_cen[0], dec_cen[0]
+        else:
+            ra_cen = pointing.ra
+            dec_cen = pointing.dec
+
         # Get WCS for the focal plane
         if wcs_fp == None:
-            wcs_fp = self.focal_plane.getTanWCS(pointing.ra, pointing.dec, pointing.img_pa, chip.pix_scale)
+            wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale)
 
         # Create chip Image
         chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
@@ -92,14 +131,18 @@ class Observation(object):
         elif chip.survey_type == "spectroscopic":
             flat_normal = np.ones_like(chip.img.array)
             if self.config["ins_effects"]["flat_fielding"] == True:
-                print("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID, flush=True)
-                print(chip.img.bounds, flush=True)
+                # print("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID, flush=True)
+                # print(chip.img.bounds, flush=True)
+                chip_output.logger.info("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID)
+                msg = str(chip.img.bounds)
+                chip_output.logger.info(msg)
                 flat_img = Effects.MakeFlatSmooth(
                     chip.img.bounds,
                     int(self.config["random_seeds"]["seed_flat"]))
                 flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array)
             if self.config["ins_effects"]["shutter_effect"] == True:
-                print("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID, flush=True)
+                # print("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID, flush=True)
+                chip_output.logger.info("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID)
                 shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735,
                                                     dt=1E-3)  # shutter effect normalized image for this chip
                 flat_normal = flat_normal*shuttimg
@@ -115,7 +158,7 @@ class Observation(object):
 
         if pointing.pointing_type == 'MS':
             # Load catalogues and templates
-            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir)
+            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, logger=chip_output.logger)
             self.nobj = len(self.cat.objs)
 
             # Loop over objects
@@ -123,8 +166,11 @@ class Observation(object):
             bright_obj = 0
             dim_obj = 0
             for j in range(self.nobj):
+                
+                # (DEBUG)
                 # if j >= 10:
                 #     break
+
                 obj = self.cat.objs[j]
                 if obj.type == 'star' and self.config["run_option"]["galaxy_only"]:
                     continue
@@ -144,7 +190,8 @@ class Observation(object):
                         norm_filt=norm_filt,
                     )
                 except Exception as e:
-                    print(e)
+                    # print(e)
+                    chip_output.logger.error(e)
                     continue
 
                 # Exclude very bright/dim objects (for now)
@@ -171,11 +218,13 @@ class Observation(object):
                         # TODO: every object with individual shear from input catalog(s)
                         obj.g1, obj.g2 = self.g1_field[j], self.g2_field[j]
                     except:
-                        print("failed to load external shear.")
+                        # print("failed to load external shear.")
+                        chip_output.logger.error("failed to load external shear.")
                         pass
                 elif self.config["shear_setting"]["shear_type"] == "catalog":
                     pass
                 else:
+                    chip_output.logger.error("Unknown shear input")
                     raise ValueError("Unknown shear input")
 
                 pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False)
@@ -224,7 +273,8 @@ class Observation(object):
                         # print("object omitted", flush=True)
                         continue
                 except Exception as e:
-                    print(e)
+                    # print(e)
+                    chip_output.logger.error(e)
                     pass
                 # Unload SED:
                 obj.unload_SED()
@@ -233,7 +283,8 @@ class Observation(object):
             del psf_model
             del self.cat
 
-        print("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
+        # print("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
+        chip_output.logger.info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
 
         # Detector Effects
         # ===========================================================
@@ -249,7 +300,8 @@ class Observation(object):
             pointing_ID=pointing.id,
             timestamp_obs=pointing.timestamp,
             pointing_type=pointing.pointing_type,
-            sky_map=sky_map, tel = self.tel)
+            sky_map=sky_map, tel = self.tel,
+            logger=chip_output.logger)
         
         if pointing.pointing_type == 'MS':
             datetime_obs = datetime.fromtimestamp(pointing.timestamp)
@@ -290,12 +342,16 @@ class Observation(object):
             hdu1 = fits.HDUList([hdu1, hdu2])
             fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits')
             hdu1.writeto(fname, output_verify='ignore', overwrite=True)
-            print("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
-            print("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
-            print("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
+            # print("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
+            # print("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
+            # print("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
+            chip_output.logger.info("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
+            chip_output.logger.info("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
+            chip_output.logger.info("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
         del chip.img
 
-        print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
+        # print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
+        chip_output.logger.info("check running:2: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
 
     def runExposure_MPI_PointingList(self, pointing_list, shear_cat_file=None, chips=None, use_mpi=False):
         if use_mpi:
@@ -334,7 +390,7 @@ class Observation(object):
 
                 chip = run_chips[ichip]
                 filt = run_filts[ichip]
-                print("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid), flush=True)
+                # print("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid), flush=True)
                 chip_output = ChipOutput(
                     config=self.config, 
                     focal_plane=self.focal_plane, 
@@ -345,6 +401,7 @@ class Observation(object):
                     pointing_ID=pointing_ID,  
                     subdir=sub_img_dir,
                     prefix=prefix)
+                chip_output.logger.info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid))
                 self.run_one_chip(
                     chip=chip, 
                     filt=filt, 
@@ -352,3 +409,6 @@ class Observation(object):
                     pointing=pointing,
                     cat_dir=self.path_dict["cat_dir"])
                 print("finished running chip#%d..."%(chip.chipID), flush=True)
+                chip_output.logger.info("finished running chip#%d..."%(chip.chipID))
+                for handler in chip_output.logger.handlers[:]:
+                    chip_output.logger.removeHandler(handler)

ObservationSim/_util.py

@@ -104,6 +104,7 @@ def make_run_dirs(work_dir, run_name, pointing_list):
                 os.makedirs(subImgDir, exist_ok=True)
             except OSError:
                 pass
+    return imgDir
 
 def imgName(tt=0):
     ut = datetime.utcnow()

config/config_NGP.yaml

@@ -12,14 +12,14 @@
 # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/"
 work_dir: "/public/home/fangyuedong/temp/CSST/workplace/"
 data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
-run_name: "NGP_20220327"
+run_name: "NGP_test"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
 # - pointing_dir: null
 # - pointing_file: null
-pointing_dir: null
-pointing_file: null
+pointing_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+pointing_file: "pointing_test_NGP_3.dat"
 
 # Whether to use MPI
 run_option:
@@ -34,7 +34,7 @@ run_option:
   out_cat_only: NO
 
   # Only simulate stars?
-  star_only: NO
+  star_only: YES
 
   # Only simulate galaxies?
   galaxy_only: NO
@@ -48,7 +48,7 @@ obs_setting:
   # "Photometric": simulate photometric chips only
   # "Spectroscopic": simulate slitless spectroscopic chips only
   # "All": simulate full focal plane
-  survey_type: "All"
+  survey_type: "Photometric"
   
   # Exposure time [seconds]
   exp_time: 150.
@@ -73,7 +73,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: null
+  run_pointings: [1]
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
@@ -92,7 +92,7 @@ obs_setting:
 # Default path settings for WIDE survey simulation
 input_path:
   cat_dir: "OnOrbitCalibration/CTargets20211231"
-  star_cat: "CT-NGP_V2.2.hdf5"
+  star_cat: "CT-NGP_r1.8_G28.hdf5"
   galaxy_cat: "galaxyCats_r_10.0_healpix_shift_192.859500_27.128300.hdf5"
 
 SED_templates_path:

run_sim.py

@@ -2,6 +2,7 @@ from ObservationSim.ObservationSim import Observation
 from ObservationSim._util import parse_args, make_run_dirs, generate_pointing_list
 import os
 import yaml
+import shutil
 
 import gc
 gc.enable()
@@ -47,7 +48,8 @@ def run_sim(Catalog):
     pointing_list = generate_pointing_list(config=config, pointing_filename=config['pointing_file'], data_dir=config['pointing_dir'])
 
     # Make the main output directories
-    make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], pointing_list=pointing_list)
+    run_dir = make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], pointing_list=pointing_list)
+    shutil.copy(args.config_file, run_dir)
 
     # Initialize the simulation
     obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir'])
@@ -63,6 +65,10 @@ if __name__=='__main__':
     # from Catalog.Catalog_example import Catalog_example
     # run_sim(Catalog=Catalog_example)
     
-    # To run cycle-3 simulation
-    from Catalog.C3Catalog import C3Catalog
-    run_sim(Catalog=C3Catalog)
+    # # To run cycle-3 simulation
+    # from Catalog.C3Catalog import C3Catalog
+    # run_sim(Catalog=C3Catalog)
+
+    # To run calibration field NGP simulation
+    from Catalog.NGPCatalog import NGPCatalog
+    run_sim(Catalog=NGPCatalog)