Merge branch 'develop' into 'master'

merge develop into master (for release v2.0)

See merge request csst_sim/csst-simulation!13

Catalog/C3Catalog.py → Catalog/C6_Catalog.py

@@ -5,15 +5,21 @@ 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 seds, sed_assign, extAv, tag_sed, getObservedSED
+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
+
 try:
     import importlib.resources as pkg_resources
 except ImportError:
@@ -22,52 +28,90 @@ except ImportError:
 
 NSIDE = 128
 
-class C3Catalog(CatalogBase):
-    def __init__(self, config, chip, pointing, **kwargs):
+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
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **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"]
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
 
-        if "logger" in kwargs:
-            self.logger = kwargs["logger"]
-        else:
-            self.logger = None
+        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))
-        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.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"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["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"])
+        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 = os.path.join(config["data_dir"], config["catalog_options"]["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.))
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
         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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
-        vertices = spherical_to_cartesian(1., dec, ra)
-        self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        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)
@@ -78,7 +122,8 @@ class C3Catalog(CatalogBase):
         if obj.type == "star":
             return self.normF_star
         elif obj.type == "galaxy" or obj.type == "quasar":
-            return self.normF_galaxy
+            # return self.normF_galaxy
+            return None
         else:
             return None
 
@@ -86,18 +131,24 @@ class C3Catalog(CatalogBase):
         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)
+        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_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_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)
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
 
         # Apply astrometric modeling
         # in C3 case only aberration
-        ra_arr = gals['ra_true'][:]
-        dec_arr = gals['dec_true'][:]
+        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()
@@ -105,7 +156,7 @@ class C3Catalog(CatalogBase):
             pmdec_list = np.zeros(ngals).tolist()
             rv_list = np.zeros(ngals).tolist()
             parallax_list = [1e-9] * ngals
-            dt = datetime.fromtimestamp(self.pointing.timestamp)
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
             date_str = dt.date().isoformat()
             time_str = dt.time().isoformat()
             ra_arr, dec_arr = on_orbit_obs_position(
@@ -122,84 +173,119 @@ class C3Catalog(CatalogBase):
                 input_vx=self.pointing.sat_vx,
                 input_vy=self.pointing.sat_vy,
                 input_vz=self.pointing.sat_vz,
-                input_epoch="J2015.5",
+                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_true'][igals]
-            param['dec_orig'] = gals['dec_true'][igals]
-            param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
-            if param['mag_use_normal'] >= 26.5:
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            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_true'][igals]
+
+            param['z'] = gals['redshift'][igals]
             param['model_tag'] = 'None'
-            param['gamma1'] = 0
-            param['gamma2'] = 0
-            param['kappa'] = 0
+            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['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            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
-            # 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:
+
+            # 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
+
             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
 
             # NOTE: this cut cannot be put before the SED type has been assigned
             if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
                 continue
 
+            # TEMP
             self.ids += 1
             # param['id'] = self.ids
-            param['id'] = gals['galaxyID'][igals]
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(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)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 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['sourceID'])
         # Apply astrometric modeling
-        # in C3 case only aberration
         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 = np.zeros(nstars).tolist()
-            pmdec_list = np.zeros(nstars).tolist()
-            rv_list = np.zeros(nstars).tolist()
-            parallax_list = [1e-9] * nstars
-            dt = datetime.fromtimestamp(self.pointing.timestamp)
+            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(
@@ -216,23 +302,30 @@ class C3Catalog(CatalogBase):
                 input_vx=self.pointing.sat_vx,
                 input_vy=self.pointing.sat_vy,
                 input_vz=self.pointing.sat_vz,
-                input_epoch="J2015.5",
+                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]
-            if param['mag_use_normal'] >= 26.5:
-                continue
+            # 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]
@@ -242,31 +335,117 @@ class C3Catalog(CatalogBase):
             param['z'] = 0.0
             param['star'] = 1   # Star
             obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            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=nAGNs,
+                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 iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
             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"]:
+        
+        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')['catalog']
             for pix in self.pix_list:
+                try:
                     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']
+                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)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
                     gals = gals_cat[str(pix)]
-                self._load_gals(gals, pix_id=pix)
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
                     del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+            try:
+                self._load_AGNs()
+            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))
-        del self.avGal
-
 
     def load_sed(self, obj, **kwargs):
         if obj.type == 'star':
@@ -278,22 +457,58 @@ class C3Catalog(CatalogBase):
                 feh=obj.param['feh']
             )
         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=self.tempSed_gal[obj.sed_type],
+                    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.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
         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)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
         y = speci(lamb)
-        # erg/s/cm2/A --> photo/s/m2/A
+        # 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)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
         del wave
         del flux
         return sed

Catalog/C6_fits_Catalog.py

+import os
+import galsim
+import random
+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, Stamp
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+import astropy.io.fits as fitsio
+from ObservationSim.MockObject._util import seds, sed_assign, extAv
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+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
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        # (TEST)
+        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"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["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["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 = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        ###mock_stamp_START
+        if "stamp_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["stamp_cat"] and config["catalog_options"]["stamp_yes"]:
+            stamp_file = config["catalog_options"]["input_path"]["stamp_cat"]
+            self.stamp_path = os.path.join(self.cat_dir, stamp_file)
+            self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir="/share/simudata/CSSOSDataProductsSims/data/Templates/Galaxy/") #only for test
+        ###mock_stamp_END
+        
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        # 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
+            return None
+        ###mock_stamp_START
+        elif obj.type == "stamp":
+            #return self.normF_galaxy  ###normalize_filter for stamp
+            return None
+        ###mock_stamp_END
+        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):
+        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_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        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]
+            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['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            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']
+
+            # 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
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 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['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.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]
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            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)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            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=nAGNs,
+                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 iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+            self.objs.append(obj)
+
+    ###mock_stamp_START
+    def _load_stamps(self, stamps, pix_id=None):
+        nstamps = len(stamps['filename'])
+        self.rng_sedGal = random.Random()
+        self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
+        self.ud = galsim.UniformDeviate(pix_id)
+
+        for istamp in range(nstamps):
+            fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
+            hdu=fitsio.open(fitsfile)
+
+            param = self.initialize_param()
+            param['id']   = hdu[0].header['index'] #istamp
+            param['star'] = 3      # Stamp type in .cat file
+            param['lensGalaxyID'] = hdu[0].header['lensGID']
+            param['ra'] = hdu[0].header['ra']
+            param['dec']= hdu[0].header['dec']
+            param['pixScale']= hdu[0].header['pixScale']
+            #param['srcGalaxyID'] = hdu[0].header['srcGID']
+            #param['mu']= hdu[0].header['mu']
+            #param['PA']= hdu[0].header['PA']
+            #param['bfrac']= hdu[0].header['bfrac']
+            #param['z']= hdu[0].header['z']
+            param['mag_use_normal'] = 22 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst']
+
+            assert(stamps['lensGID'][istamp] == param['lensGalaxyID'])
+
+            # Apply astrometric modeling
+            # in C3 case only aberration
+            param['ra_orig'] = param['ra']
+            param['dec_orig']= param['dec']
+            if self.config["obs_setting"]["enable_astrometric_model"]:
+                ra_list = [param['ra']] #ra_arr.tolist()
+                dec_list= [param['dec']] #dec_arr.tolist()
+                pmra_list = np.zeros(1).tolist()
+                pmdec_list = np.zeros(1).tolist()
+                rv_list = np.zeros(1).tolist()
+                parallax_list = [1e-9] * 1
+                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=1,
+                    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
+                )
+                param['ra'] = ra_arr[0]
+                param['dec']= dec_arr[0]
+
+            # 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'] = 0.0
+            param['redden'] = 0
+
+            #param["CSSTmag"]= True
+            #param["mag_r"] = 20.
+            #param['']
+            ###more keywords for stamp###
+            param['image'] = hdu[0].data
+            param['image'] = param['image']/(np.sum(param['image']))
+            obj = Stamp(param)
+            self.objs.append(obj)
+    ###mock_stamp_END
+
+    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')['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)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+            try:
+                self._load_AGNs()
+            except Exception as e:
+                traceback.print_exc()
+                self.logger.error(str(e))
+                print(e)
+
+        ###mock_stamp_START
+        if "stamp_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["stamp_cat"] and self.config["catalog_options"]["stamp_yes"]:
+            stamps_cat = h5.File(self.stamp_path, 'r')['Stamps']
+            for pix in self.pix_list:
+                try:
+                    stamps = stamps_cat[str(pix)]
+                    self._load_stamps(stamps, pix_id=pix)
+                    del stamps
+                except Exception as e:
+                    self.logger.error(str(e))
+                    print(e)
+        ###mock_stamp_END
+
+        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']
+            )
+        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.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        ###mock_stamp_START
+        elif obj.type == 'stamp':
+            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]
+        ###mock_stamp_END
+        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)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/Catalog_example.py

@@ -6,7 +6,7 @@ from scipy import interpolate
 
 from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
 
-class Catalog_example(CatalogBase):
+class Catalog(CatalogBase):
     """An user customizable class for reading in catalog(s) of objects and SEDs.
     
     NOTE: must inherit the "CatalogBase" abstract class
@@ -32,8 +32,7 @@ class Catalog_example(CatalogBase):
     load_norm_filt(obj):
         load the filter throughput for the input catalog's photometric system.
     """
-
-    def __init__(self, config, chip, pointing, **kwargs):
+    def __init__(self, config, chip, **kwargs):
         """Constructor method.
         
         Parameters
@@ -41,7 +40,7 @@ class Catalog_example(CatalogBase):
         config : dict
             configuration dictionary which is parsed from the input YAML file
         chip: ObservationSim.Instrument.Chip
-            a ObservationSim.Instrument.Chip instance, can be used to identify the band etc.
+            an ObservationSim.Instrument.Chip instance, can be used to identify the band etc.
         **kwargs : dict
             other needed input parameters (in key-value pairs), please modify corresponding
             initialization call in "ObservationSim.py" as you need.
@@ -52,12 +51,11 @@ class Catalog_example(CatalogBase):
         """
         
         super().__init__()
-        self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
         self.chip = chip
-        self.pointing = pointing
-        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
-            star_file = config["input_path"]["star_cat"]
-            star_SED_file = config["SED_templates_path"]["star_SED"]
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"]
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["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)
         # NOTE: must call _load() method here to read in all objects
@@ -100,20 +98,23 @@ class Catalog_example(CatalogBase):
                 NOTE: if that filter is not the corresponding CSST filter, the 
                 load_norm_filt(obj) function must be implemented to load the filter
                 throughput of that particular photometric system
-            param["theta"] : float
-                the position angle (in degrees)
             param["bfrac"] : float
                 the bulge fraction
-            param["hlr_bulge] : float
+            param["hlr_bulge"] : float
                 the half-light-radius of the bulge
-            param["hlr_disk] : float
+            param["hlr_disk"] : float
                 the half-light-radius of the disk
-             param["ell_bulge] : float
-                the ellipticity of the bulge
-             param["ell_disk] : float
-                the ellipticity of the disk
-             param["ell_tot] : float
-                the total ellipticity
+            param["e1_bulge"], param["e2_bulge"] : float
+                the ellipticity of the bulge components
+            param["e1_disk"], param["e2_disk"] : float
+                the ellipticity of the disk components
+            (Optional parameters):
+            param['disk_sersic_idx']: float
+                Sersic index for galaxy disk component
+            param['bulge_sersic_idx']: float
+                Sersic index for galaxy bulge component
+            param['g1'], param['g2']: float
+                Reduced weak lensing shear components (valid for shear type: catalog)
         the model of ObservationSim.MockObject.Galaxy class requires:
             Currently a Quasar is modeled as a point source, just like a Star.
 
@@ -122,6 +123,9 @@ class Catalog_example(CatalogBase):
 
         NOTE: All constructed objects should be appened to "self.objs".
 
+        NOTE: Any other parameters can also be set within "param" dict:
+            Used to calculate required quantities and/or SEDs etc.
+
         Parameters
         ----------
         **kwargs : dict

Catalog/FGS_Catalog.py

+import os
+import galsim
+import random
+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
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+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
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        # (TEST)
+        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"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["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["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 = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        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_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_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        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]
+            # 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['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            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
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 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['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.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]
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            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)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            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=nAGNs,
+                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 iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    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')['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)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+            try:
+                self._load_AGNs()
+            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']
+            )
+        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.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        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':
+        if obj.type == 'galaxy' or 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)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/NGPCatalog.py

@@ -22,11 +22,11 @@ except ImportError:
 
 NSIDE = 128
 
-class NGPCatalog(CatalogBase):
+class Catalog(CatalogBase):
     def __init__(self, config, chip, pointing, chip_output, **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"]
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
 
         self.chip_output = chip_output
         self.logger = chip_output.logger
@@ -40,19 +40,19 @@ class NGPCatalog(CatalogBase):
         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"]
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["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"]
+        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_file = config["catalog_options"]["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.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["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.))
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
         else:
             self.rotation = 0.
 
@@ -111,7 +111,7 @@ class NGPCatalog(CatalogBase):
             pmdec_list = np.zeros(ngals).tolist()
             rv_list = np.zeros(ngals).tolist()
             parallax_list = [1e-9] * ngals
-            dt = datetime.fromtimestamp(self.pointing.timestamp)
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
             date_str = dt.date().isoformat()
             time_str = dt.time().isoformat()
             ra_arr, dec_arr = on_orbit_obs_position(
@@ -140,12 +140,12 @@ class NGPCatalog(CatalogBase):
             param['ra_orig'] = gals['ra_true'][igals]
             param['dec_orig'] = gals['dec_true'][igals]
             param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
-            if param['mag_use_normal'] >= 26.5:
-                continue
+            # 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['g1'] = 0
+            param['g2'] = 0
             param['kappa'] = 0
             param['delta_ra'] = 0
             param['delta_dec'] = 0
@@ -167,8 +167,13 @@ class NGPCatalog(CatalogBase):
 
             # 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)]
+            # param['redden'] = self.tempRed_gal[param['sed_type']]
+            # param['av'] = self.avGal[int(self.ud()*self.nav)]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
             if param['sed_type'] <= 5:
                 param['av'] = 0.0
                 param['redden'] = 0
@@ -211,7 +216,7 @@ class NGPCatalog(CatalogBase):
             pmdec_list = pmdec_arr.tolist()
             rv_list = rv_arr.tolist()
             parallax_list = parallax_arr.tolist()
-            dt = datetime.fromtimestamp(self.pointing.timestamp)
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
             date_str = dt.date().isoformat()
             time_str = dt.time().isoformat()
             ra_arr, dec_arr = on_orbit_obs_position(
@@ -245,8 +250,8 @@ class NGPCatalog(CatalogBase):
             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
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
             self.ids += 1
             # param['id'] = self.ids
             param['id'] = stars['sourceID'][istars]
@@ -269,7 +274,7 @@ class NGPCatalog(CatalogBase):
         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"]:
+        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')['catalog']
             for pix in self.pix_list:
                 try:
@@ -279,7 +284,7 @@ class NGPCatalog(CatalogBase):
                 except Exception as e:
                     self.logger.error(str(e))
                     print(e)
-        if "galaxy_cat" in self.config["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
+        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"]:
             gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
             for pix in self.pix_list:
                 try:
@@ -316,8 +321,6 @@ class NGPCatalog(CatalogBase):
         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)
         lamb = np.arange(2000, 18001 + 0.5, 0.5)
         y = speci(lamb)
         # erg/s/cm2/A --> photo/s/m2/A

Catalog/testCat_fits.py

+import os
+import galsim
+import random
+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, Stamp
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+import astropy.io.fits as fitsio
+from ObservationSim.MockObject._util import seds, sed_assign, extAv
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+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 Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        # (TEST)
+        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 "stamp_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["stamp_cat"] and config["catalog_options"]["stamp_yes"]:
+            stamp_file = config["catalog_options"]["input_path"]["stamp_cat"]
+            self.stamp_path = os.path.join(self.cat_dir, stamp_file)
+            #self.stamp_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["stamp_SED"]) ###shoule be stamp-SED
+            #self._load_SED_lib_stamps() ###shoule be stamp-SED
+            self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir="/share/simudata/CSSOSDataProductsSims/data/Templates/Galaxy/") #only for test
+        
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        # 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 == "stamp":
+            #return self.normF_galaxy  ###normalize_filter for stamp
+            return None
+        else:
+            return None
+
+    def _load_stamps(self, stamps, pix_id=None):
+        nstamps = len(stamps['filename'])
+        self.rng_sedGal = random.Random()
+        self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
+        self.ud = galsim.UniformDeviate(pix_id)
+
+        for istamp in range(nstamps):
+            fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
+            hdu=fitsio.open(fitsfile)
+
+            param = self.initialize_param()
+            param['id']   = hdu[0].header['index'] #istamp
+            param['star'] = 3      # Stamp type in .cat file
+            param['ra'] = hdu[0].header['ra']
+            param['dec']= hdu[0].header['dec']
+            param['pixScale']= hdu[0].header['pixScale']
+            #param['srcGalaxyID'] = hdu[0].header['srcGID']
+            #param['mu']= hdu[0].header['mu']
+            #param['PA']= hdu[0].header['PA']
+            #param['bfrac']= hdu[0].header['bfrac']
+            #param['z']= hdu[0].header['z']
+            param['mag_use_normal'] = hdu[0].header['mag_g'] #gals['mag_true_g_lsst']
+
+            # Apply astrometric modeling
+            # in C3 case only aberration
+            param['ra_orig'] = param['ra']
+            param['dec_orig']= param['dec']
+            if self.config["obs_setting"]["enable_astrometric_model"]:
+                ra_list = [param['ra']] #ra_arr.tolist()
+                dec_list= [param['dec']] #dec_arr.tolist()
+                pmra_list = np.zeros(1).tolist()
+                pmdec_list = np.zeros(1).tolist()
+                rv_list = np.zeros(1).tolist()
+                parallax_list = [1e-9] * 1
+                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=1,
+                    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
+                )
+                param['ra'] = ra_arr[0]
+                param['dec']= dec_arr[0]
+
+            # 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'] = 0.0
+            param['redden'] = 0
+
+            #param["CSSTmag"]= True
+            #param["mag_r"] = 20.
+            #param['']
+            ###more keywords for stamp###
+            param['image'] = hdu[0].data
+            param['image'] = param['image']/(np.sum(param['image']))
+            obj = Stamp(param)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+
+        if "stamp_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["stamp_cat"] and self.config["catalog_options"]["stamp_yes"]:
+            stamps_cat = h5.File(self.stamp_path, 'r')['Stamps']
+            for pix in self.pix_list:
+                try:
+                    stamps = stamps_cat[str(pix)]
+                    self._load_stamps(stamps, pix_id=pix)
+                    del stamps
+                except Exception as e:
+                    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 == 'stamp':
+            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(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'))
+        
+        del wave
+        del flux
+        return sed

Catalog/testCat_galaxy.py

+import os
+import galsim
+import random
+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, Stamp
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+import astropy.io.fits as fitsio
+from ObservationSim.MockObject._util import seds, sed_assign, extAv
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+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
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        # (TEST)
+        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 "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and config["catalog_options"]["galaxy_yes"]:
+            self.galaxy_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["galaxy_cat"])
+
+            self.galaxy_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and config["catalog_options"]["galaxy_yes"]:
+            self.AGN_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["AGN_cat"])
+
+            self.AGN_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        # 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 == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    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_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        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
+            # if igals < 3447:
+            #     continue
+            
+            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]
+            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
+
+            # if param['mag_use_normal'] >= 26.5:
+            #     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['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            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']
+
+            # 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
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 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_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            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=nAGNs,
+                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 iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    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 "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and self.config["catalog_options"]["galaxy_yes"]:
+            # for i in range(76):
+            #     file_path = os.path.join(self.galaxy_path, "galaxies%04d_C6.h5"%(i))
+            #     gals_cat = h5.File(file_path, 'r')['galaxies']
+            #     for pix in self.pix_list:
+            #         try:
+            #             gals = gals_cat[str(pix)]
+            #             self._load_gals(gals, pix_id=pix, cat_id=i)
+            #             del gals
+            #         except Exception as e:
+            #             traceback.print_exc()
+            #             self.logger.error(str(e))
+            #             print(e)
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    # # (TEST C6):
+                    # if pix != 35421 or bundleID != 523:
+                    #     continue
+                    gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+        #if "AGN_cat" in self.config["input_path"] and self.config["input_path"]["AGN_cat"] and not self.config["run_option"]["star_only"]:
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and self.config["catalog_options"]["galaxy_yes"]:
+            try:
+                self._load_AGNs()
+            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))
+
+    # (TEST)
+    # def convert_mags(self, obj):
+    #     spec = np.matmul(obj.coeff, self.pcs.T)
+    #     lamb = self.lamb_gal * U.angstrom
+    #     unit_sed = ((lamb * lamb)/constants.c*U.erg/U.second/U.cm**2/U.angstrom).to(U.jansky)
+    #     unit_sed = unit_sed.value
+    #     lamb = lamb.value
+    #     lamb *= (1 + obj.z)
+    #     spec *= (1 + obj.z)
+    #     mags = -2.5 * np.log10(unit_sed*spec) + 8.9
+    #     mags += self.cosmo.distmod(obj.z).value
+    #     return lamb, mags
+
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'galaxy' or obj.type == 'quasar':
+            # dist_L_pc = (1 + obj.z) * comoving_dist(z=obj.z)[0]
+            # factor = (10 / dist_L_pc)**2
+            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.SED_AGN[int(obj.id)] * factor * 1e-17
+                flux = self.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        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)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/testCat_star.py

+import os
+import galsim
+import random
+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, Stamp
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+import astropy.io.fits as fitsio
+from ObservationSim.MockObject._util import seds, sed_assign, extAv
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+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 Catalog(CatalogBase):
+    """An user customizable class for reading in catalog(s) of objects and SEDs.
+    
+    NOTE: must inherit the "CatalogBase" abstract class
+
+    ...
+    
+    Attributes
+    ----------
+    cat_dir : str
+        a directory that contains the catalog file(s)
+    star_path : str
+        path to the star catalog file
+    star_SED_path : str
+        path to the star SED data
+    objs : list
+        a list of ObservationSim.MockObject (Star, Galaxy, or Quasar)
+        NOTE: must have "obj" list when implement your own Catalog
+    
+    Methods
+    ----------
+    load_sed(obj, **kwargs):
+        load the corresponding SED data for one object
+    load_norm_filt(obj):
+        load the filter throughput for the input catalog's photometric system.
+    """
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        """Constructor method.
+        
+        Parameters
+        ----------
+        config : dict
+            configuration dictionary which is parsed from the input YAML file
+        chip: ObservationSim.Instrument.Chip
+            an ObservationSim.Instrument.Chip instance, can be used to identify the band etc.
+        pointing: ObservationSim.Config.Pointing
+            an ObservationSim.Config.Pointing instance, can be used to configure the astrometry module
+        chip_output: ObservationSim.Config.ChipOutput
+            an ObservationSim.Config.ChipOutput instance, can be used to setup the output format
+        filt: ObservationSim.Instrument.Filter
+            an ObservationSim.Instrument.Filter instance, can be used to identify the filter type
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+        
+        Returns
+        ----------
+        None
+        """
+        
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        #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 config["catalog_options"]["star_yes"]:
+            self.star_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["star_cat"])
+
+            self.star_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["star_SED"])
+            self._load_SED_lib_star()
+        
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        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 = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %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_ouptut_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]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        # 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):
+        """Load the corresponding thourghput for the input magnitude "param["mag_use_normal"]".
+
+        NOTE: if the input magnitude is already in CSST magnitude, simply return None
+
+        Parameters
+        ----------
+        obj : ObservationSim.MockObject
+            the object to get thourghput data for
+
+        Returns
+        ----------
+        norm_filt : Astropy.Table
+            the throughput Table with two columns (namely, "WAVELENGTH", "SENSITIVITY"):
+                norm_filt["WAVELENGTH"] : wavelengthes in Angstroms
+                norm_filt["SENSITIVITY"] : efficiencies
+        """
+        if obj.type == "star":
+            return self.normF_star
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    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.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]
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            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)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], 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):
+        """Read in all objects in from the catalog file(s).
+
+        This is a must implemented method which is used to read in all objects, and
+        then convert them to ObservationSim.MockObject (Star, Galaxy, or Quasar).
+        
+        Currently, 
+        the model of ObservationSim.MockObject.Star class requires:
+            param["star"] : int
+                specify the object type: 0: galaxy, 1: star, 2: quasar
+            param["id"] : int
+                ID number of the object
+            param["ra"] : float
+                Right ascension (in degrees)
+            param["dec"] : float
+                Declination (in degrees)
+            param["mag_use_normal"]: float
+                the absolute magnitude in a particular filter
+                NOTE: if that filter is not the corresponding CSST filter, the 
+                load_norm_filt(obj) function must be implemented to load the filter
+                throughput of that particular photometric system
+        
+        the model of ObservationSim.MockObject.Galaxy class requires:
+            param["star"] : int
+                specify the object type: 0: galaxy, 1: star, 2: quasar
+            param["id"] : int
+                ID number of the object
+            param["ra"] : float
+                Right ascension (in degrees)
+            param["dec"] : float
+                Declination (in degrees)
+            param["mag_use_normal"]: float
+                the absolute magnitude in a particular filter
+                NOTE: if that filter is not the corresponding CSST filter, the 
+                load_norm_filt(obj) function must be implemented to load the filter
+                throughput of that particular photometric system
+            param["bfrac"] : float
+                the bulge fraction
+            param["hlr_bulge"] : float
+                the half-light-radius of the bulge
+            param["hlr_disk"] : float
+                the half-light-radius of the disk
+            param["e1_bulge"], param["e2_bulge"] : float
+                the ellipticity of the bulge components
+            param["e1_disk"], param["e2_disk"] : float
+                the ellipticity of the disk components
+            (Optional parameters):
+            param['disk_sersic_idx']: float
+                Sersic index for galaxy disk component
+            param['bulge_sersic_idx']: float
+                Sersic index for galaxy bulge component
+            param['g1'], param['g2']: float
+                Reduced weak lensing shear components (valid for shear type: catalog)
+        the model of ObservationSim.MockObject.Galaxy class requires:
+            Currently a Quasar is modeled as a point source, just like a Star.
+
+        NOTE: To construct an object, according to its type, just call:
+            Star(param), Galaxy(param), or Quasar(param)
+
+        NOTE: All constructed objects should be appened to "self.objs".
+
+        NOTE: Any other parameters can also be set within "param" dict:
+            Used to calculate required quantities and/or SEDs etc.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+
+        Returns
+        ----------
+        None
+        """
+        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"]:
+        if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and self.config["catalog_options"]["star_yes"]:
+            star_cat = h5.File(self.star_path, 'r')['stars']
+            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 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))
+
+    def load_sed(self, obj, **kwargs):
+        """Load the corresponding SED data for a particular obj.
+
+        Parameters
+        ----------
+        obj : ObservationSim.MockObject
+            the object to get SED data for
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+
+        Returns
+        ----------
+        sed : Astropy.Table
+            the SED Table with two columns (namely, "WAVELENGTH", "FLUX"):
+                sed["WAVELENGTH"] : wavelength in Angstroms
+                sed["FLUX"] : fluxes in photons/s/m^2/A
+            NOTE: the range of wavelengthes must at least cover [2450 - 11000] Angstorms
+        """
+
+        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']
+            )
+        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'))
+        
+        del wave
+        del flux
+        return sed

ObservationSim/Config/ChipOutput.py

@@ -18,43 +18,39 @@ class ChipOutput(object):
         else:
             self.ra_cen = config["obs_setting"]["ra_center"]
             self.dec_cen = config["obs_setting"]["dec_center"]
-        exp_name = imgKey0 + "_%s_%s.fits"
-        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_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat"
+        self.chipLabel = focal_plane.getChipLabel(chip.chipID)
+        self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
 
         self.subdir = subdir
 
         # Setup logger for each chip
-        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log"
+        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
         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)
         logging.getLogger('numba').setLevel(logging.WARNING)
-        # formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
         formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(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 "
-        hdr4 = "pm_ra pm_dec RV parallax"
-
-        fmt1  = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s "
-        fmt2  = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f "
-        fmt3 = "%2d %8.4f %8.4f "
-        fmt4 = "%15.8f %15.8f %15.8f %15.8f"
-
-        self.hdr = hdr1 + hdr2 + hdr3 + hdr4
-        self.fmt = fmt1 + fmt2 + fmt3 + fmt4
+        hdr2 = "pm_ra pm_dec RV parallax"
+        fmt1  = "%20s %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s "
+        fmt2 = "%15.8f %15.8f %15.8f %15.8f"
+        self.hdr = hdr1 + hdr2
+        self.fmt = fmt1 + fmt2
 
         self.logger.info("pointing_type = %s\n"%(pointing_type))
     
+    def Log_info(self, message):
+        print(message)
+        self.logger.info(message)
+    
+    def Log_error(self, message):
+        print(message)
+        self.logger.error(message)
     
     def update_ouptut_header(self, additional_column_names=""):
         self.hdr += additional_column_names
@@ -70,16 +66,14 @@ class ChipOutput(object):
     def cat_add_obj(self, obj, pos_img, pos_shear):
         # ximg = pos_img.x - self.chip.bound.xmin + 1.0
         # yimg = pos_img.y - self.chip.bound.ymin + 1.0
-        # print('-------------debug-----------------')
-        # print('from global',ximg, yimg)
+        # self.logger.info('-------------debug-----------------')
+        # self.logger.info('from global',ximg, yimg)
         ximg = obj.real_pos.x + 1.0
         yimg = obj.real_pos.y + 1.0
-        # print('from loacl',ximg, yimg)
+        # self.logger.info('from loacl',ximg, yimg)
         
         line = self.fmt%(
             obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, 
-            obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2,
-            obj.sed_type, obj.av, obj.redden,
             obj.pmra, obj.pmdec, obj.rv, obj.parallax)
         line += obj.additional_output_str
         if not line.endswith("\n"):

ObservationSim/Config/Config.py

@@ -17,14 +17,9 @@ def config_dir(config, work_dir=None, data_dir=None):
         path_dict["data_dir"] =os.path.join(path_dict["work_dir"], "data/")
     else:
         path_dict["data_dir"] = data_dir
-    # Data sub-catalogs
-    # Object catalog direcotry
-    path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], config["input_path"]["cat_dir"])
     # PSF data directory
     if config["psf_setting"]["psf_model"] == "Interp":
         path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"])
-    if config["ins_effects"]["field_dist"] == True:
-        path_dict["fd_path"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["fd_path"])
 
     return path_dict
 

ObservationSim/Config/Header/ImageHeader.py

@@ -13,7 +13,10 @@ import random
 import os
 import sys
 import astropy.coordinates as coord
+from astropy.coordinates import SkyCoord
+from astropy.wcs.utils import fit_wcs_from_points
 from astropy.time import Time
+from astropy import wcs
 
 def chara2digit(char):
     """ Function to judge and convert characters to digitals
@@ -164,7 +167,8 @@ def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232
 
 
 ##9232 9216  898 534 1309 60 -40  -23.4333
-def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'):
+def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra_ref = 60, dec_ref = -40, pa = -23.433, pixel_scale = 0.074, pixel_size=1e-2, 
+    rotate_chip=0., filter = 'GI', row_num = None, col_num = None, xcen = None, ycen = None):
 
     """ Creat a wcs frame for CCST with multiple extensions
 
@@ -172,26 +176,30 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
     ----------
 
     """
+    r_dat = OrderedDict()
+    r_dat['EQUINOX'] = 2000.0
+    r_dat['WCSDIM'] = 2.0
+    r_dat['CTYPE1'] = 'RA---TAN'
+    r_dat['CTYPE2'] = 'DEC--TAN'
+    r_dat['CRVAL1'] = ra_ref
+    r_dat['CRVAL2'] = dec_ref
     
     flag_x = [0, 1, -1, 1, -1]
     flag_y = [0, 1, 1, -1, -1]
     flag_ext_x = [0,-1,1,-1,1]
     flag_ext_y = [0,-1,-1,1,1]
+    pa_aper = pa
+
+    if (row_num is not None) and (col_num is not None):
+
         x_num = 6
         y_num = 5
         detector_num = x_num*y_num
 
-
         detector_size_x = xlen
         detector_size_y = ylen
         gap_y = gapy
         gap_x = [gapx1,gapx2]
-    ra_ref = ra
-    dec_ref = dec
-
-    pa_aper = pa
-
-    pixel_size = psize
 
         gap_x1_num = 3
         gap_x2_num = 2
@@ -202,37 +210,18 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
 
         gap_x_map = np.array([[0,0,0,0,0],[gap_x[0],gap_x[1],gap_x[1],gap_x[1],gap_x[1]],[gap_x[1],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[1]],[gap_x[1],gap_x[1],gap_x[1],gap_x[1],gap_x[0]]])
         
-
-    # frame_array = np.empty((5,6),dtype=np.float64)
-    # print(x_center,y_center)
-    
         j = row_num
         i = col_num
         # ccdnum = str((j-1)*5+i)
 
         x_ref, y_ref = detector_size_x*i + sum(gap_x_map[0:i,j-1]) - detector_size_x/2. , (detector_size_y+gap_y)*j-gap_y-detector_size_y/2
         
-    # print(i,j,x_ref,y_ref,ra_ref,dec_ref)
-
-    r_dat = OrderedDict()
-
-    # name = []
-    # value = []
-    # description = []
-    
         for k in range(1,2):
             
-        cd = np.array([[ pixel_size,  0], [0, pixel_size]])/3600.*flag_x[k]
+            cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[k]
             cd_rot = rotate_CD_matrix(cd, pa_aper)
 
             # f = open("CCD"+ccdnum.rjust(2,'0')+"_extension"+str(k)+"_wcs.param","w")
-
-        r_dat['EQUINOX'] = 2000.0
-        r_dat['WCSDIM'] = 2.0
-        r_dat['CTYPE1'] = 'RA---TAN'
-        r_dat['CTYPE2'] = 'DEC--TAN'
-        r_dat['CRVAL1'] = ra_ref
-        r_dat['CRVAL2'] = dec_ref
             r_dat['CRPIX1'] = flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center)
             r_dat['CRPIX2'] = flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center)
             r_dat['CD1_1'] = cd_rot[0,0]
@@ -241,7 +230,6 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
             r_dat['CD2_2'] = cd_rot[1,1]
 
             if filter in ['GU', 'GV', 'GI']:
-            from astropy import wcs
 
                 w = wcs.WCS(naxis=2)
                 w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
@@ -271,13 +259,10 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
                         y_pixs[i_pix] = y
 
                         pix_coor = np.array([x, y])
-                    sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w)
+                        sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, xlen=xlen, ylen=ylen, w=w)
                         # print(sc1[0,0],sc1[0,1])
                         sky_coors.append((sc1[0, 0], sc1[0, 1]))
 
-            from astropy.coordinates import SkyCoord
-            from astropy.wcs.utils import fit_wcs_from_points
-
                 wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
                                             world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN')
 
@@ -296,17 +281,67 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
                 r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
                 r_dat['CRVAL2'] = wcs_new.wcs.crval[1]
 
-        return r_dat
+    elif (xcen is not None) and (ycen is not None):
+        xcen, ycen = xcen/pixel_size, ycen/pixel_size
+        x1, y1 = xcen - xlen/2., ycen - ylen/2.
+        r_dat['CRPIX1'] = -x1
+        r_dat['CRPIX2'] = -y1
+
+        cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
+        cd_rot = rotate_CD_matrix(cd, pa_aper)
+        r_dat['CD1_1'] = cd_rot[0,0]
+        r_dat['CD1_2'] = cd_rot[0,1]
+        r_dat['CD2_1'] = cd_rot[1,0]
+        r_dat['CD2_2'] = cd_rot[1,1]
+
+        w = wcs.WCS(naxis=2)
+        w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
+        w.wcs.cd = cd_rot
+        w.wcs.crval = [ra_ref, dec_ref]
+        w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]
 
+        sn_x = 30
+        sn_y = 30
+        x_pixs = np.zeros(sn_y * sn_x)
+        y_pixs = np.zeros(sn_y * sn_x)
+        xpixs_line = np.linspace(1, xlen, sn_x)
+        ypixs_line = np.linspace(1, ylen, sn_y)
+        sky_coors = []
+        for n1, y in enumerate(ypixs_line):
+            for n2, x in enumerate(xpixs_line):
+                i_pix = n1 * sn_x + n2
+                x_pixs[i_pix] = x
+                y_pixs[i_pix] = y
+                pix_coor = np.array([x, y])
+                sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=rotate_chip, xlen=xlen, ylen=ylen, w=w)
+                sky_coors.append((sc1[0, 0], sc1[0, 1]))
+        wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]), 
+                                      world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), 
+                                      projection='TAN')
+        r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
+        r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
+        r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
+        r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
+        r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
+        r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]
+
+        r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
+        r_dat['CRVAL2'] = wcs_new.wcs.crval[1]
+
+    else:
+        raise ValueError('In function WCS_def(): Either (row_num, col_num) or (xcen, ycen, pixel_size) should be given')
+    
+    return r_dat
 
 
 
-def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.]):
+#TODO project_cycle is temporary， is not in header defined, delete in future
+def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, chip_name="01"):
 
     # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
 
     
-    k = (row_num-1)*6+col_num
+    # k = (row_num-1)*6+col_num
     # ccdnum = str(k)
 
     g_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/global_header.header'
@@ -325,48 +360,64 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     h_prim = fits.Header()
     h_prim = fits.Header.fromfile(g_header_fn)
 
-    h_prim['PIXSIZE1'] = xlen
-    h_prim['PIXSIZE2'] = ylen
-
-    h_prim['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6]
-    h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
-    h_prim['DATE-OBS'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6]
-    h_prim['TIME-OBS'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
-    h_prim['DETECTOR'] = 'CCD'+CCDID[k-1].rjust(2,'0')
-    h_prim['RA_OBJ'] = ra
-    h_prim['DEC_OBJ'] = dec
-    h_prim['OBJECT'] = '1'+ pointNum.rjust(8,'0')
-    h_prim['OBSID'] = '1'+ pointNum.rjust(8,'0')
-    h_prim['TELFOCUS'] = 'f/14'
+    # h_prim['PIXSIZE1'] = xlen
+    # h_prim['PIXSIZE2'] = ylen
+
+    h_prim['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
+    # h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
+    h_prim['DATE-OBS'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
+    # h_prim['TIME-OBS'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
+    # h_prim['DETECTOR'] = 'CHIP'+CCDID[k-1].rjust(2,'0')
+    h_prim['OBJ_RA'] = ra
+    h_prim['OBJ_DEC'] = dec
+
+    obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
+
+    OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
+
+    h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
+    h_prim['OBSID'] = OBS_id
+    # h_prim['TELFOCUS'] = 'f/14'
     h_prim['EXPTIME'] = exptime
     
     # Define file types
-    file_type = {'SCI':'sci', 'BIAS':'zero', 'DARK':'dark', 'FLAT':'flat', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray'}
+    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'}
     h_prim['FILETYPE'] = file_type[im_type]
 
     co = coord.SkyCoord(ra, dec, unit='deg')
 
-    ra_hms = format(co.ra.hms.h, '02.0f') + ':' + format(co.ra.hms.m, '02.0f') + ':' + format(co.ra.hms.s, '05.2f')
-    dec_hms = format(co.dec.dms.d, '02.0f') + ':' + format(abs(co.dec.dms.m), '02.0f') + ':' + format(abs(co.dec.dms.s),
-                                                                                                      '05.2f')
+    ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m, '02.0f') + format(co.ra.hms.s, '02.0f')
+    dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m), '02.0f') + format(abs(co.dec.dms.s), '02.0f')
+    h_prim['TARGET'] = ra_hms + '+' + dec_hms
+    #
+    # h_prim['RA_NOM'] = ra_hms
+    # h_prim['DEC_NOM'] = dec_hms
 
-    h_prim['RA_NOM'] = ra_hms
-    h_prim['DEC_NOM'] = dec_hms
-    h_prim['PIXSCAL1'] = psize
-    h_prim['PIXSCAL2'] = psize
+    h_prim['RA_PNT0'] = ra
+    h_prim['DEC_PNT0'] = dec
+    h_prim['RA_PNT1'] = ra
+    h_prim['DEC_PNT1'] = dec
 
-    ttt = h_prim['DATE'] + 'T' + h_prim['TIME']
+
+
+    # h_prim['PIXSCAL1'] = pixel_scale
+    # h_prim['PIXSCAL2'] = pixel_scale
+
+    ttt = h_prim['DATE']
     tstart = Time(ttt)
     h_prim['EXPSTART'] = round(tstart.mjd, 5)
+    h_prim['CABSTART'] = h_prim['EXPSTART']
     # tend = Time(tstart.cxcsec + h_prim['EXPTIME'], format="cxcsec")
     tend = Time(tstart.mjd + h_prim['EXPTIME']/86400., format="mjd")
     h_prim['EXPEND'] = round(tend.mjd, 5)
+    h_prim['CABEND'] = h_prim['EXPEND']
 
     file_start_time = '20' + date[0:6] + time_obs[0:6]
     end_time_str = str(tend.datetime)
     file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
-    h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_1' + pointNum.rjust(8, '0') + '_' + CCDID[
-        k - 1].rjust(2, '0') + '_L0_1'
+    # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
+    #     k - 1].rjust(2, '0') + '_L0_V01'
+    h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
 
 
     h_prim['POSI0_X'] = sat_pos[0]
@@ -376,17 +427,18 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     h_prim['VELO0_X'] = sat_vel[0]
     h_prim['VELO0_Y'] = sat_vel[1]
     h_prim['VELO0_Z'] = sat_vel[2]
-    h_prim['RA_PNT0'] = ra_hms
-    h_prim['DEC_PNT0'] = dec_hms
+    # h_prim['RA_PNT0'] = ra_hms
+    # h_prim['DEC_PNT0'] = dec_hms
     
     # Get version of CSSTSim Package
     from pkg_resources import get_distribution
-    h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
-
+    # h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
+    h_prim['FITSCREA'] = get_distribution("CSSTSim").version
 
     return h_prim
 
-def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1, extName='SCI'):
+def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2, 
+                            extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None):
 
     e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'
     f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
@@ -401,21 +453,58 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2
     s = s.strip("\n")
     CCDID = s.split()
 
-    k = (row_num - 1) * 6 + col_num
+    # k = (row_num - 1) * 6 + col_num
 
     h_ext = fits.Header.fromfile(e_header_fn)
 
-    h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
-    h_ext['FILTER'] = filters[k-1]
+    # h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
+    # h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
+    # h_ext['FILTER'] = filters[k-1]
+    h_ext['CCDCHIP'] = chip.chip_name
+    h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0')
+    h_ext['FILTER'] = chip.filter_type
     h_ext['NAXIS1'] = xlen
     h_ext['NAXIS2'] = ylen
     h_ext['EXTNAME'] = extName
     h_ext['GAIN1'] = gain
-    h_ext['RDNOISE1'] = readout
-    h_ext['CCDCHIP'] =  'ccd' + CCDID[k-1].rjust(2,'0')
-    h_ext['POS_ANG'] = pa
-    header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize,
-                         row_num=row_num, col_num=col_num, filter = h_ext['FILTER'])
+    h_ext['GAIN2'] = gain
+    h_ext['GAIN3'] = gain
+    h_ext['GAIN4'] = gain
+    h_ext['GAIN5'] = gain
+    h_ext['GAIN6'] = gain
+    h_ext['GAIN7'] = gain
+    h_ext['GAIN8'] = gain
+    h_ext['GAIN9'] = gain
+    h_ext['GAIN10'] = gain
+    h_ext['GAIN11'] = gain
+    h_ext['GAIN12'] = gain
+    h_ext['GAIN13'] = gain
+    h_ext['GAIN14'] = gain
+    h_ext['GAIN15'] = gain
+    h_ext['GAIN16'] = gain
+    h_ext['RDNOIS1'] = readout
+    h_ext['RDNOIS2'] = readout
+    h_ext['RDNOIS3'] = readout
+    h_ext['RDNOIS4'] = readout
+    h_ext['RDNOIS5'] = readout
+    h_ext['RDNOIS6'] = readout
+    h_ext['RDNOIS7'] = readout
+    h_ext['RDNOIS8'] = readout
+    h_ext['RDNOIS9'] = readout
+    h_ext['RDNOIS10'] = readout
+    h_ext['RDNOIS11'] = readout
+    h_ext['RDNOIS12'] = readout
+    h_ext['RDNOIS13'] = readout
+    h_ext['RDNOIS14'] = readout
+    h_ext['RDNOIS15'] = readout
+    h_ext['RDNOIS16'] = readout
+
+    h_ext['PIXSCAL1'] = pixel_scale
+    h_ext['PIXSCAL2'] = pixel_scale
+
+    # 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, 
+                         rotate_chip=chip.rotate_angle, filter = h_ext['FILTER'], row_num=row_num, col_num=col_num, xcen = xcen, ycen = ycen)
 
     h_ext['CRPIX1'] = header_wcs['CRPIX1']
     h_ext['CRPIX2'] = header_wcs['CRPIX2']

ObservationSim/Config/Header/extension_header.header

-XTENSION= 'IMAGE   '           / Image extension                                BITPIX  =                   16 / array data type                                NAXIS   =                    2 / number of array dimensions                     NAXIS1  =                 9216 / length of first array axis                     NAXIS2  =                 9232 / length of second array axis                    PCOUNT  =                    0                                                  GCOUNT  =                    1                                                  EXTNAME = 'SCI     '                                                            EXTVER  =                    1                                                  BSCALE  =                    1                                                  BZERO   =                32768                                                  BUNIT   = 'ADU     '                                                            COMMENT ==================================================================      COMMENT CCD chip information                                                    COMMENT ==================================================================      NCHAN   =                   16 / Number of readout channels                     NCHAN1  =                    8 / Number of horizintal channels                  NCHAN2  =                    2 / Number of vertical channels                    DETSIZE = '[1:9216:9232]'      / detector size                                  BIASSEC = '[9216:9217,9232:9234]' / bias section                                CCDCHIP =                    6 / CCD chip ID                                    CCDLABEL= 'y-1     '           / CCD chip label                                 PSCAN1  =                   27 / horizontal prescan width                       PSCAN2  =                    8 / vertical prescan height                        OSCAN1  =                   16 / horizontal overscan width                      OSCAN2  =                    8 / vertical overscan height                       FILTER  = 'y       '           / filter name                                    COMMENT ==================================================================      COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS                          COMMENT ==================================================================      WCSDIM  =                    2 / Number of World Coordinate System axes         EQUINOX =               2000.0 / Epoch (year)                                   CRPIX1  =             -10017.0 / Coordinate reference pixel of x                CRPIX2  =              24876.0 / Coordinate reference pixel of y                CRVAL1  =            62.228226 / Coordinate reference value of x                CRVAL2  =           -42.316932 / Coordinate reference value of y                CTYPE1  = 'RA---TAN'           / the coordinate type                            CTYPE2  = 'DEC--TAN'           / the coordinate type                            CD1_1   = 1.88602083707394E-05 / partial of first axis coordinate of x          CD2_1   = -8.1745583617600E-06 / partial of first axis coordinate of y          CD1_2   = 8.17455836176000E-06 / partial of second axis coordinate of x         CD2_2   = 1.88602083707394E-05 / partial of second axis coordinate of y         POS_ANG =                 30.0 / Positon angle of detector array                COMMENT ==================================================================      COMMENT Readout information                                                     COMMENT ==================================================================      GAIN1   =                  1.1 / CCD gain                                       RDNOISE1=                  5.0 / read noise                                     READTIME=                 40.0 / read time (sec)                                RDSPEED =                 10.0 / read speed (in MHz)                            CHIPTEMP=               -100.0 / chip temperature (in K)                        COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= '''abcde'''          / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
+XTENSION= 'IMAGE   '           / extension type                                 BITPIX  =                   16 / bits per data value                            NAXIS   =                    2 / number of data axes                            NAXIS1  =                 9216 / length of first array axis                     NAXIS2  =                 9232 / length of second array axis                    PCOUNT  =                    0                                                  GCOUNT  =                    1                                                  EXTNAME = 'SCI     '                                                            EXTVER  =                    1                                                  BSCALE  =                    1                                                  BZERO   =                32768                                                  BUNIT   = 'ADU     '            / physical unit of array values                 COMMENT ==================================================================      COMMENT Detector information                                                    COMMENT ==================================================================      DETECTOR= 'CCD'                 / detector name                                 DETTEMP1=                   173 / detector temperature at EXPSTART              DETTEMP2=                   173 / detector temperature at EXPEND                DETTEMP3=                   173 / detector temperature when readout is finished DETSIZE = '9560x9264'           / detector size                                 DATASEC = '9216x9232'           / data size                                     PIXSCAL1=                 0.074 / pixel scale for axis 1                        PIXSCAL2=                 0.074 / pixel scale for axis 2                        PIXSIZE1=                    10 / pixel size  for axis 1                        PIXSIZE2=                    10 / pixel size  for axis 1                        COMMENT ==================================================================      COMMENT CCD chip information                                                    COMMENT ==================================================================      CCDCHIP =                     8 / CCD chip ID                                   CCDLABEL= 'y-1'                 / CCD chip label                                FILTER  = 'y'                   / filter name                                   NCHAN   =                    16 / number of readout channels                    NCHAN1  =                     8 / number of horizintal channels                 NCHAN2  =                     2 / number of vertical channels                   PSCAN1  =                    27 / horizontal prescan width, per readout channel PSCAN2  =                     8 / horizontal prescan width, per readout channel OSCAN1  =                    16 / horizontal overscan width,per readout channel OSCAN2  =                    16 / horizontal overscan width,per readout channel COMMENT ==================================================================      COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS                          COMMENT ==================================================================      WCSAXES =                    2 / number of World Coordinate System axes         CRPIX1  =             -10017.0 / x-coordinate of reference pixel                CRPIX2  =              24876.0 / y-coordinate of reference pixel                CRVAL1  =            62.228226 / first axis value at reference pixel            CRVAL2  =           -42.316932 / second axis value at reference pixel           CTYPE1  = 'RA---TAN'           / the coordinate type for the first axis         CTYPE2  = 'DEC--TAN'           / the coordinate type for the second axis        CD1_1   = 1.88602083707394E-05 / partial of first axis coordinate w.r.t. x      CD2_1   = -8.1745583617600E-06 / partial of first axis coordinate w.r.t. x      CD1_2   = 8.17455836176000E-06 / partial of second axis coordinate w.r.t. x     CD2_2   = 1.88602083707394E-05 / partial of second axis coordinate w.r.t. x     OTHERS  =                 ''   /                                                COMMENT ==================================================================      COMMENT Readout information                                                     COMMENT ==================================================================      GAIN1   =                  1.1 / CCD gain (channel 1)                           GAIN2   =                  1.1 / CCD gain (channel 2)                           GAIN3   =                  1.1 / CCD gain (channel 3)                           GAIN4   =                  1.1 / CCD gain (channel 4)                           GAIN5   =                  1.1 / CCD gain (channel 5)                           GAIN6   =                  1.1 / CCD gain (channel 6)                           GAIN7   =                  1.1 / CCD gain (channel 7)                           GAIN8   =                  1.1 / CCD gain (channel 8)                           GAIN9   =                  1.1 / CCD gain (channel 9)                           GAIN10  =                  1.1 / CCD gain (channel 10)                          GAIN11  =                  1.1 / CCD gain (channel 11)                          GAIN12  =                  1.1 / CCD gain (channel 12)                          GAIN13  =                  1.1 / CCD gain (channel 13)                          GAIN14  =                  1.1 / CCD gain (channel 14)                          GAIN15  =                  1.1 / CCD gain (channel 15)                          GAIN16  =                  1.1 / CCD gain (channel 16)                          RDNOIS1 =                  5.0 / read noise  (channel 1)                        RDNOIS2 =                  5.0 / read noise  (channel 2)                        RDNOIS3 =                  5.0 / read noise  (channel 3)                        RDNOIS4 =                  5.0 / read noise  (channel 4)                        RDNOIS5 =                  5.0 / read noise  (channel 5)                        RDNOIS6 =                  5.0 / read noise  (channel 6)                        RDNOIS7 =                  5.0 / read noise  (channel 7)                        RDNOIS8 =                  5.0 / read noise  (channel 8)                        RDNOIS9 =                  5.0 / read noise  (channel 9)                        RDNOIS10=                  5.0 / read noise  (channel 10)                       RDNOIS11=                  5.0 / read noise  (channel 11)                       RDNOIS12=                  5.0 / read noise  (channel 12)                       RDNOIS13=                  5.0 / read noise  (channel 13)                       RDNOIS14=                  5.0 / read noise  (channel 14)                       RDNOIS15=                  5.0 / read noise  (channel 15)                       RDNOIS16=                  5.0 / read noise  (channel 16)                       RDSPEED =                 10.0 / read speed (in MHz)                            COMMENT ==================================================================      COMMENT Shutter information                                                     COMMENT ==================================================================      SHUTSTAT=                    T / shutter status                                 SHTOPEN0=                  0.0 / shutter open time (begin)                      SHTOPEN1=                  0.0 / shutter open time (end)                        SHTCLOS0=                  0.0 / shutter close time (begin)                     SHTCLOS1=                  0.0 / shutter close time (end)                       COMMENT ==================================================================      COMMENT LED information                                                         COMMENT ==================================================================      LEDSTAT1= '0000000000000000000000000000' / main LED status                      LEDEXP1 =                  0.0 / main LED flash time (s)                        LEDSTAT2= '0000000000000000000000000000' / backup LED status                    LEDEXP2 =                  0.0 / backup LED flash time (s)                      LEDTEMPA=                173.0 / LED temperature (main LED in K)                LEDTEMPB=                173.0 / LED temperature (backup LED in K)              COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= '''abcde'''          / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          

ObservationSim/Config/Header/global_header.header

-SIMPLE  =                    T / conforms to FITS standard                      BITPIX  =                    8 / array data type                                NAXIS   =                    0 / number of array dimensions                     EXTEND  =                    T                                                  NEXTEND =                    1                                                  GROUPS  =                    F                                                  DATE    = '2021-03-04'         / Date this file was written                     TIME    = '09:30:00'           / Time this file was written                     FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / Name of file                 FILETYPE= 'sci     '           / type of data found in data data file           TELESCOP= 'CSST    '           / telescope used to acquire data                 INSTRUME= 'MSC     '           / identifier for instrument used to acquire data RADECSYS= 'ICRS    '           / frame of reference of coordinates              EQUINOX =               2000.0 / Default equinox                                COMMENT ==================================================================      COMMENT Object information                                                      COMMENT ==================================================================      OBJECT  = '00000000'           / name of the object observed                    TARGET  = '+000000000000'      / Name of the target (RA & Dec of the object)    OBSID   = '00000000'           / ID of this observation                         RA_OBJ  =            62.228226 / R.A. of the object (degrees)                   DEC_OBJ =           -42.316932 / declination of the object (degrees)            COMMENT ==================================================================      COMMENT Exposure information                                                    COMMENT ==================================================================      SUNANGLE=                 30.0 / angle between sun and direction perpendicular  MOONANGL=                 30.0 / angle between moon and direction perpendicular SUN_ALT =                 30.0 / Sun altitude                                   REFFRAME= 'CSSTGSC-1.0'        / guide star catalog version                     DATE-OBS= '2021-03-04'         / UTC date of start of observation (yyyy-mm-dd)  TIME-OBS= '09:30:00'           / UTC time of start of observation (hh:mm:ss)    EXPSTART=              58849.0 / exposure start time (MJD)                      EXPEND  =              58849.0 / exposure end time (MJD)                        EXPTIME =                150.0 / exposure duration                              EPOCH   =               2021.0 / coordinate epoch                               COMMENT ==================================================================      COMMENT Detector information                                                    COMMENT ==================================================================      DETECTOR= 'CCD06   '           / deterctor name                                 PIXSCAL1=                0.074 / Pixel scale for axis 1                         PIXSCAL2=                0.074 / Pixel scale for axis 2                         PIXSIZE1=               9216.0 / Pixel size for axis 1                          PIXSIZE2=               9232.0 / Pixel size for aixs 2                          CCDTEMP =               -100.0 / CCD chip / dewar temperature (in K)            COMMENT ==================================================================      COMMENT Telescope information                                                   COMMENT ==================================================================      TELSTAT =                    T / Tracking status                                TELTEMP1=               -100.0 / Telescope temperature (in K)                   POSI0_X =                  0.0 / The orbital position in X (shuttr open)        POSI0_Y =                  0.0 / The orbital position in Y (shuttr open)        POSI0_Z =                  0.0 / The orbital position in Z (shuttr open)        VELO0_X =                  0.0 / The orbital velocity in X (shuttr open)        VELO0_Y =                  0.0 / The orbital velocity in Y (shuttr open)        VELO0_Z =                  0.0 / The orbital velocity in Z (shuttr open)        RA_PNT0 = '00:00:00.00'        / R.A. at shutter open time                      DEC_PNT0= '00:00:00.00'        / declination at shutter open time               POSI1_X =                  0.0 / The orbital position in X (shuttr close)       POSI1_Y =                  0.0 / The orbital position in Y (shuttr close)       POSI1_Z =                  0.0 / The orbital position in Z (shuttr close)       VELO1_X =                  0.0 / The orbital velocity in X (shuttr close)       VELO1_Y =                  0.0 / The orbital velocity in Y (shuttr close)       VELO1_Z =                  0.0 / The orbital velocity in Z (shuttr close)       RA_PNT1 =                  0.0 / R.A. at shutter close time                     DEC_PNT1=                  0.0 / declination at shutter close time              RA_NOM  = '00:00:00.00'        / nominal proposed R.A. of the pointing          DEC_NOM = '00:00:00.00'        / nominal proposed DEC of the pointing           TELFOCUS= 'abcdefg '           / telescope focus                                AOSTAT  =                    T / Active opticas status                          HOODSTAT=                    T / Lens hood status                               COMMENT ==================================================================      COMMENT Shutter information                                                     COMMENT ==================================================================      SHUTHWV = 'shutter-h-1.0'      / shutter hardware version                       SHUTSWV = 'shutter-s-1.0'      / shutter software version                       SHUTSTAT=                    T / shutter status                                 SHTOPEN0=                  0.0 / shutter open time (begin)                      SHTOPEN1=                  0.0 / shutter open time (end)                        SHTCLOS0=                  0.0 / shutter close time (begin)                     SHTCLOS1=                  0.0 / shutter close time (end)                       COMMENT ==================================================================      COMMENT Guide information                                                       COMMENT ==================================================================      STARNO1 =                  1.0 / guide star number                              STARPOX1=                  0.0 / guide star average position                    STARPOY1=                  0.0 / guide star average position                    SECMOM1 =                  0.0 / guide star second momentum                     GAFLAG1 =                    T / guide flag indicate whether adjust position    STARNO2 =                  1.0 / guide star number                              STARPOX2=                  0.0 / guide star average position                    STARPOY2=                  0.0 / guide star average position                    SECMOM2 =                  0.0 / guide star second momentum                     GAFLAG2 =                    T / guide flag indicate whether adjust position    STARNO3 =                  1.0 / guide star number                              STARPOX3=                  0.0 / guide star average position                    STARPOY3=                  0.0 / guide star average position                    SECMOM3 =                  0.0 / guide star second momentum                     GAFLAG3 =                    T / guide flag indicate whether adjust position    STARNO4 =                  1.0 / guide star number                              STARPOX4=                  0.0 / guide star average position                    STARPOY4=                  0.0 / guide star average position                    SECMOM4 =                  0.0 / guide star second momentum                     GAFLAG4 =                    T / guide flag indicate whether adjust position    COMMENT ==================================================================      COMMENT LED information                                                         COMMENT ==================================================================      LEDSTAT1= '0000000000000000000000000000' / LED status                           LEDEXP1 =                  0.0 / LED flash time (second)                        LEDSTAT2= '0000000000000000000000000000' / LED status                           LEDEXP2 =                  0.0 / LED flash time (second)                        LEDTEMPA=               -100.0 / LED temperature (main LED in K)                LEDTEMPB=               -100.0 / LED temperature (backup LED in K)              COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= 'abcdefg '           / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
+SIMPLE  =                    T / conforms to FITS standard                      BITPIX  =                    8 / array data type                                NAXIS   =                    0 / number of array dimensions                     EXTEND  =                    T                                                  NEXTEND =                    1                                                  GROUPS  =                    F                                                  DATE    = '2021-03-04T09:30:00'/ date this file was written                     FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name                    FILETYPE= 'SCI     '           / observation type                               TELESCOP= 'CSST    '           / telescope used to acquire data                 INSTRUME= 'MSC     '           / instrument used to acquire data                RADECSYS= 'ICRS    '           / frame of reference of coordinates              EQUINOX =               2000.0 /                                                FITSCREA= 'C6'                 / FITS create software version                   COMMENT ==================================================================      COMMENT Object information                                                      COMMENT ==================================================================      OBJECT  = '00000000'           / object name                                    TARGET  = '+000000000000'      / target name (hhmmss+ddmmss)                    OBSID   = '00000000'           / observation ID                                 OBJ_RA  =            62.228226 / R.A. of the object (degrees)                   OBJ_DEC =           -42.316932 / declination of the object (degrees)            COMMENT ==================================================================      COMMENT Telescope information                                                   COMMENT ==================================================================      REFFRAME= 'CSSTGSC-1.0'        / guide star catalog version                     DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss)  SATESWV = '0001'               / software version in the satellite              EXPSTART=              59130.5 / exposure start time (MJD)                      CABSTART=              59130.5 / (MJD)                                          SUNANGL0=                 50.0 / angle between sun and optical axis at CABSTART MOONANG0=                 30.0 / angle moon and opt axis at CABSTART at CABST   TEL_ALT0=                 20.0 / angle opt axis and the ground-piston at CABST  POS_ANG0=                 20.0 / angle bwt y axis and the North Pole at CABST   POSI0_X =                  0.0 / the orbital position in X at CABSTART          POSI0_Y =                  0.0 / the orbital position in Y at CABSTART          POSI0_Z =                  0.0 / the orbital position in Z at CABSTART          VELO0_X =                  0.0 / the orbital velocity in X at CABSTART          VELO0_Y =                  0.0 / the orbital velocity in Y at CABSTART          VELO0_Z =                  0.0 / the orbital velocity in Z at CABSTART          EULER0_1=                  0.0 / Euler angle 1 at  CABSTART                     EULER0_2=                  0.0 / Euler angle 2 at  CABSTART                     EULER0_3=                  0.0 / Euler angle 3 at  CABSTART                     RA_PNT0 =                  0.0 / RA of the pointing (degrees) at CABSTART       DEC_PNT0=                  0.0 / DEC of the pointing (degrees) at CABSTART      EXPEND  =                  0.0 / exposure end time (MJD)                        CABEND  =                  0.0 / (MJD)                                          SUNANGL1=                 50.0 / angle between sun and optical axis at CABEND   MOONANG1=                 30.0 / angle btw moon and optical axis at CABEND      TEL_ALT1=                 20.0 / angle opt axis and the ground-piston at CAEND  POS_ANG1=                 20.0 / angle bwt y axis and the North Pole at CABEND  POSI1_X =                  0.0 / the orbital position in X at CABEND            POSI1_Y =                  0.0 / the orbital position in Y at CABEND            POSI1_Z =                  0.0 / the orbital position in Z at CABEND            VELO1_X =                  0.0 / The orbital velocity in X at CABEND            VELO1_Y =                  0.0 / The orbital velocity in Y at CABEND            VELO1_Z =                  0.0 / The orbital velocity in Z at CABEND            EULER1_1=                  0.0 / Euler angle 1 at  CABEND                       EULER1_2=                  0.0 / Euler angle 2 at  CABEND                       EULER1_3=                  0.0 / Euler angle 3 at  CABEND                       RA_PNT1 =                  0.0 / RA of the pointing (degrees) at CABEND         DEC_PNT1=                  0.0 / DEC of the pointing (degrees) at CABEND        EXPTIME =                150.0 / exposure duration                              EPOCH   =                 2000 / coordinate epoch                               COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= 'abcdefg '           / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             

ObservationSim/Instrument/Chip/Chip.py

 import galsim
 import os
 import numpy as np
+import pickle
+import json
 from astropy.table import Table
 from numpy.random import Generator, PCG64
 from astropy.io import fits
@@ -9,6 +11,7 @@ from datetime import datetime
 from ObservationSim.Instrument.Chip import Effects as effects
 from ObservationSim.Instrument.FocalPlane import FocalPlane
 from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
+from ObservationSim.Instrument._util import rotate_conterclockwise
 
 try:
     import importlib.resources as pkg_resources
@@ -21,30 +24,65 @@ class Chip(FocalPlane):
         # Get focal plane (instance of paraent class) info
         # TODO: use chipID to config individual chip?
         super().__init__()
-        self.npix_x = 9216
-        self.npix_y = 9232
-        self.read_noise = 5.0   # e/pix
-        self.dark_noise = 0.02  # e/pix/s
-        self.pix_scale  = 0.074 # pixel scale
-        self.gain = float(config["ins_effects"]["gain"])
-        self.bias_level = float(config["ins_effects"]["bias_level"])
-        self.overscan   = 1000
-        self.exptime    = float(config["obs_setting"]["exp_time"])   # second
-        self.dark_exptime = float(config["ins_effects"]['dark_exptime'])
-        self.flat_exptime = float(config["ins_effects"]['flat_exptime'])
-        self.readout_time = float(config["ins_effects"]['readout_time'])
-        self.full_well = int(config["ins_effects"]["full_well"])
+        # self.npix_x = 9216
+        # self.npix_y = 9232
+        # self.pix_scale  = 0.074 # pixel scale
+
+        self._set_attributes_from_config(config)
 
         self.logger = logger
 
         # A chip ID must be assigned
         self.chipID = int(chipID)
-        self._getChipRowCol()
+        self.chip_name = str(chipID).rjust(2, '0')
 
         # Get corresponding filter info
         self.filter_id, self.filter_type = self.getChipFilter()
         self.survey_type = self._getSurveyType()
 
+        # [TODO]
+        if self.filter_type != "FGS":
+            self._getChipRowCol()
+
+        # Set the relavent specs for FGS detectors
+        # [TODO]
+        try:
+            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition:
+                with open(chip_definition, "r") as f:
+                    chip_dict = json.load(f)[str(self.chipID)]
+        except AttributeError:
+            with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition:
+                with open(chip_definition, "r") as f:
+                    chip_dict = json.load(f)[str(self.chipID)]
+        for key in chip_dict:
+            setattr(self, key, chip_dict[key])
+        if self.filter_type == "FGS":
+            if config["ins_effects"]["field_dist"] == False:
+                self.fdModel = None
+            else:
+                fgs_name = self.chip_name[0:4]
+                try:
+                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
+                        with open(field_distortion, "rb") as f:
+                            self.fdModel = pickle.load(f)
+                except AttributeError:
+                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
+                        with open(field_distortion, "rb") as f:
+                            self.fdModel = pickle.load(f)
+        else:
+            # Get the corresponding field distortion model
+            if config["ins_effects"]["field_dist"] == False:
+                self.fdModel = None
+            else:
+                try:
+                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                        with open(field_distortion, "rb") as f:
+                            self.fdModel = pickle.load(f)
+                except AttributeError:
+                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                        with open(field_distortion, "rb") as f:
+                            self.fdModel = pickle.load(f)
+
         # Get boundary (in pix)
         self.bound = self.getChipLim()
         self.ccdEffCurve_dir = ccdEffCurve_dir
@@ -76,17 +114,25 @@ class Chip(FocalPlane):
         self.effCurve = self._getChipEffCurve(self.filter_type)
         self._getCRdata()
 
-        # Define the sensor
+        # Define the sensor model
         if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
-            self.sensor = galsim.SiliconSensor(strength=config["ins_effects"]["df_strength"], treering_func=treering_func)
+            self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
         else:
             self.sensor = galsim.Sensor()
 
         self.flat_cube = None # for spectroscopic flat field cube simulation
 
-    # def _getChipRowCol(self):
-    #     self.rowID = (self.chipID - 1) // 5 + 1
-    #     self.colID = (self.chipID - 1) % 5 + 1
+    def _set_attributes_from_config(self, config):
+        # Default setting
+        self.read_noise = 5.0   # e/pix
+        self.dark_noise = 0.02  # e/pix/s
+        self.rotate_angle = 0.
+        self.overscan   = 1000
+        # Override default values
+        for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]:
+            if key in config["ins_effects"]:
+                setattr(self, key, config["ins_effects"][key])
+
     def _getChipRowCol(self):
         self.rowID, self.colID = self.getChipRowCol(self.chipID)
 
@@ -98,20 +144,21 @@ class Chip(FocalPlane):
     def _getSurveyType(self):
         if self.filter_type in ["GI", "GV", "GU"]:
             return "spectroscopic"
-        else:
+        elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
             return "photometric"
+        # elif self.filter_type in ["FGS"]:
+        #     return "FGS"
 
     def _getChipEffCurve(self, filter_type):
         # CCD efficiency curves
         if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
-        if filter_type in ['g', 'r', 'GV']: filename = 'Astro_MB.txt'
+        if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
         if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
         # Mirror efficiency:
         # if filter_type == 'nuv': mirror_eff = 0.54
         # if filter_type == 'u': mirror_eff = 0.68
         # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
         # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
-        
         # path = os.path.join(self.ccdEffCurve_dir, filename)
         # table = Table.read(path, format='ascii')
         try:
@@ -120,14 +167,11 @@ class Chip(FocalPlane):
         except AttributeError:
             with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
                 table = Table.read(ccd_path, format='ascii')
-        # throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
         throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
         bandpass = galsim.Bandpass(throughput, wave_type='nm')
         return bandpass
 
     def _getCRdata(self):
-        # path = os.path.join(self.CRdata_dir, 'wfc-cr-attachpixel.dat')
-        # self.attachedSizes = np.loadtxt(path)
         try:
             with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path:
                 self.attachedSizes = np.loadtxt(cr_path)
@@ -138,26 +182,14 @@ class Chip(FocalPlane):
     def getChipFilter(self, chipID=None, filter_layout=None):
         """Return the filter index and type for a given chip #(chipID)
         """
-        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI"]
-        # TODO: maybe a more elegent way other than hard coded?
-        # e.g. use something like a nested dict:
+        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
         if filter_layout is not None:
             return filter_layout[chipID][0], filter_layout[chipID][1]
         if chipID == None:
             chipID = self.chipID
 
         # updated configurations
-        # if chipID>30 or chipID<1: raise ValueError("!!! Chip ID: [1,30]")
-        # if chipID in [10, 15, 16, 21]: filter_type = 'y'
-        # if chipID in [11, 20]:         filter_type = "z"
-        # if chipID in [9, 22]:           filter_type = "i"
-        # if chipID in [12, 19]:         filter_type = "u"
-        # if chipID in [7, 24]:         filter_type = "r"
-        # if chipID in [14, 13, 18, 17]:    filter_type = "nuv"
-        # if chipID in [8, 23]:         filter_type = "g"
-        # if chipID in [6, 5, 25, 26]:    filter_type = "GI"
-        # if chipID in [27, 30, 1, 4]:    filter_type = "GV"
-        # if chipID in [28, 29, 2, 3]:    filter_type = "GU"
+        if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]")
         if chipID in [6, 15, 16, 25]:  filter_type = "y"
         if chipID in [11, 20]:         filter_type = "z"
         if chipID in [7, 24]:          filter_type = "i"
@@ -168,39 +200,22 @@ class Chip(FocalPlane):
         if chipID in [1, 10, 21, 30]:  filter_type = "GI"
         if chipID in [2, 5, 26, 29]:   filter_type = "GV"
         if chipID in [3, 4, 27, 28]:   filter_type = "GU"
+        if chipID in range(31, 43):    filter_type = 'FGS'
         filter_id = filter_type_list.index(filter_type)
 
         return filter_id, filter_type
 
     def getChipLim(self, chipID=None):
         """Calculate the edges in pixel for a given CCD chip on the focal plane
-        NOTE: There are 5*4 CCD chips in the focus plane for photometric observation.
+        NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
         Parameters:
             chipID:         int
                             the index of the chip
         Returns:
             A galsim BoundsD object
         """
-        # if chipID == None:
-        #     chipID = self.chipID
-        
-        # gx = self.npix_gap_x
-        # gy1, gy2 = self.npix_gap_y
-
-        # # xlim of a given ccd chip
-        # xrem = (chipID-1)%self.nchip_x - self.nchip_x // 2
-        # xcen = (self.npix_x + gx) * xrem
-        # nx0 = xcen - self.npix_x//2 + 1
-        # nx1 = xcen + self.npix_x//2
-
-        # # ylim of a given ccd chip
-        # yrem = 2*((chipID-1)//self.nchip_x) - (self.nchip_y-1)
-        # ycen = (self.npix_y//2 + gy1//2) * yrem
-        # if chipID <= 6: ycen = (self.npix_y//2 + gy1//2) * yrem - (gy2-gy1)
-        # if chipID >= 25: ycen = (self.npix_y//2 + gy1//2) * yrem + (gy2-gy1)
-        # ny0 = ycen - self.npix_y//2 + 1
-        # ny1 = ycen + self.npix_y//2
-
+        if ((chipID is not None) and (int(chipID) <= 30)) or (self.chipID <= 30):
+            # [TODO]
             if chipID == None:
                 chipID = self.chipID
                 rowID, colID = self.rowID, self.colID
@@ -224,11 +239,25 @@ class Chip(FocalPlane):
             ycen = (self.npix_y + gy) * yrem
             ny0 = ycen - self.npix_y//2 + 1
             ny1 = ycen + self.npix_y//2
-
             return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1)
+        else:
+            xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10
+            xcen = self.x_cen / self.pix_size
+            ycen = self.y_cen / self.pix_size
+            sign_x = [-1., 1., -1., 1.]
+            sign_y = [-1., -1., 1., 1.]
+            for i in range(4):
+                x = xcen + sign_x[i] * self.npix_x / 2.
+                y = ycen + sign_y[i] * self.npix_y / 2.
+                x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
+                xmin, xmax = min(xmin, x), max(xmax, x)
+                ymin, ymax = min(ymin, y), max(ymax, y)
+            return galsim.BoundsD(xmin, xmax, ymin, ymax)
+
 
 
     def getSkyCoverage(self, wcs):
+        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
         return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)
 
 
@@ -239,14 +268,19 @@ class Chip(FocalPlane):
         bound = self.getSkyCoverage(wcs)
         return galsim.BoundsD(bound.xmin - margin, bound.xmax + margin, bound.ymin - margin, bound.ymax + margin)
 
-    def isContainObj(self, ra_obj, dec_obj, wcs=None, margin=1):
+    def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
         # magin in number of pix
+        if (ra_obj is not None) and (dec_obj is not None):
             if wcs is None:
                 wcs = self.img.wcs
             pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
+            x_image, y_image = pos_obj.x, pos_obj.y
+        elif (x_image is None) or (y_image is None):
+            raise ValueError("Either (ra_obj, dec_obj) or (x_image, y_image) should be given")
+
         xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
         ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
-        if (pos_obj.x - xmin) * (pos_obj.x - xmax) > 0.0 or (pos_obj.y - ymin) * (pos_obj.y - ymax) > 0.0:
+        if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
             return False
         return True
 
@@ -277,15 +311,15 @@ class Chip(FocalPlane):
             pointNum = str(pointing_ID),
             ra=ra_cen, 
             dec=dec_cen, 
-            psize=self.pix_scale, 
-            row_num=self.rowID, 
-            col_num=self.colID,
+            pixel_scale=self.pix_scale,
             date=date_obs,
             time_obs=time_obs,
             im_type = im_type,
-            exptime=exptime
+            exptime=exptime,
+            chip_name=str(self.chipID).rjust(2, '0')
             )
         h_ext = generateExtensionHeader(
+            chip=self,
             xlen=self.npix_x, 
             ylen=self.npix_y, 
             ra=ra_cen, 
@@ -295,9 +329,10 @@ class Chip(FocalPlane):
             readout=self.read_noise, 
             dark=self.dark_noise, 
             saturation=90000, 
-            psize=self.pix_scale, 
-            row_num=self.rowID, 
-            col_num=self.colID,
+            pixel_scale=self.pix_scale, 
+            pixel_size=self.pix_size,
+            xcen=self.x_cen,
+            ycen=self.y_cen,
             extName='raw')
         return h_prim, h_ext
 
@@ -312,7 +347,9 @@ class Chip(FocalPlane):
             time_obs=time_obs,
             exptime=exptime)
         hdu1 = fits.PrimaryHDU(header=h_prim)
+        hdu1.add_checksum()
         hdu2 = fits.ImageHDU(img.array, header=h_ext)
+        hdu2.add_checksum()
         hdu1 = fits.HDUList([hdu1, hdu2])
         fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
         hdu1.writeto(fname, output_verify='ignore', overwrite=True)
@@ -324,7 +361,7 @@ class Chip(FocalPlane):
         SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"])
         SeedDefective = int(config["random_seeds"]["seed_defective"])
         SeedCosmicRay = int(config["random_seeds"]["seed_CR"])
-        fullwell = int(config["ins_effects"]["full_well"])
+        fullwell = int(self.full_well)
         if config["ins_effects"]["add_hotpixels"] == True:
             BoolHotPix = True
         else:
@@ -342,7 +379,7 @@ class Chip(FocalPlane):
 
         # Add sky background
         if sky_map is None:
-            sky_map = filt.getSkyNoise(exptime=self.exptime)
+            sky_map = filt.getSkyNoise(exptime=exptime)
             sky_map = sky_map * np.ones_like(img.array)
             sky_map = galsim.Image(array=sky_map)
             # Apply Poisson noise to the sky map
@@ -354,7 +391,7 @@ class Chip(FocalPlane):
         elif img.array.shape != sky_map.shape:
             raise ValueError("The shape img and sky_map must be equal.")
         elif tel is not None: # If sky_map is given in flux
-            sky_map = sky_map * tel.pupil_area * self.exptime
+            sky_map = sky_map * tel.pupil_area * exptime
         if config["ins_effects"]["add_back"] == True:
             img += sky_map
         del sky_map
@@ -407,8 +444,8 @@ class Chip(FocalPlane):
                 print("  Adding Cosmic-Ray", flush=True)
             cr_map, cr_event_num = effects.produceCR_Map(
                 xLen=self.npix_x, yLen=self.npix_y, 
-                exTime=self.exptime+0.5*self.readout_time, 
-                cr_pixelRatio=0.003*(self.exptime+0.5*self.readout_time)/600.,
+                exTime=exptime+0.5*self.readout_time, 
+                cr_pixelRatio=0.003*(exptime+0.5*self.readout_time)/600.,
                 gain=self.gain, 
                 attachedSizes=self.attachedSizes,
                 seed=SeedCosmicRay+pointing_ID*30+self.chipID)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
@@ -419,7 +456,7 @@ class Chip(FocalPlane):
             del cr_map
             # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
             # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
-            datetime_obs = datetime.fromtimestamp(timestamp_obs)
+            datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
             date_obs = datetime_obs.strftime("%y%m%d")
             time_obs = datetime_obs.strftime("%H%M%S")
             self.outputCal(
@@ -432,7 +469,7 @@ class Chip(FocalPlane):
                 date_obs=date_obs,
                 time_obs=time_obs,
                 output_dir=chip_output.subdir,
-                exptime=self.exptime)
+                exptime=exptime)
             del crmap_gsimg
 
         # Apply PRNU effect and output PRNU flat file:
@@ -454,7 +491,7 @@ class Chip(FocalPlane):
 
         # Add dark current
         if config["ins_effects"]["add_dark"] == True:
-            dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(self.exptime+0.5*self.readout_time)))
+            dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
             img.addNoise(dark_noise)
 
         # Add Hot Pixels or/and Dead Pixels
@@ -474,7 +511,7 @@ class Chip(FocalPlane):
                 print("  Adding Bias level and 16-channel non-uniformity")
             if config["ins_effects"]["bias_16channel"] == True:
                 img = effects.AddBiasNonUniform16(img, 
-                    bias_level=float(config["ins_effects"]["bias_level"]), 
+                    bias_level=float(self.bias_level), 
                     nsecy = 2, nsecx=8, 
                     seed=SeedBiasNonuni+self.chipID,
                     logger=self.logger)
@@ -539,11 +576,11 @@ class Chip(FocalPlane):
                 self.logger.info("  Output N frame Bias files")
             else:
                 print("  Output N frame Bias files", flush=True)
-            NBias = int(config["ins_effects"]["NBias"])
+            NBias = int(config["output_setting"]["NBias"])
             for i in range(NBias):
                 BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                     self.npix_x, self.npix_y, 
-                    bias_level=float(config["ins_effects"]["bias_level"]), 
+                    bias_level=float(self.bias_level), 
                     ncombine=1, read_noise=self.read_noise, gain=1,
                     seed=SeedBiasNonuni+self.chipID,
                     logger=self.logger)
@@ -571,7 +608,7 @@ class Chip(FocalPlane):
 
                 # 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, logger=self.logger) + float(config["ins_effects"]["bias_level"])-5
+                    BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
 
                 BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
@@ -585,7 +622,7 @@ class Chip(FocalPlane):
                 BiasCombImg.quantize()
                 BiasCombImg = galsim.ImageUS(BiasCombImg)
                 # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
-                datetime_obs = datetime.fromtimestamp(timestamp_obs)
+                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60
@@ -608,7 +645,7 @@ class Chip(FocalPlane):
                 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"])
+            NFlat = int(config["output_setting"]["NFlat"])
             if config["ins_effects"]["add_bias"] == True:
                 biaslevel = self.bias_level
                 overscan = biaslevel-2
@@ -689,7 +726,7 @@ class Chip(FocalPlane):
                 FlatCombImg.quantize()
                 FlatCombImg = galsim.ImageUS(FlatCombImg)
                 # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
-                datetime_obs = datetime.fromtimestamp(timestamp_obs)
+                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60
@@ -717,7 +754,7 @@ class Chip(FocalPlane):
                 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"])
+            NDark = int(config["output_setting"]["NDark"])
             if config["ins_effects"]["add_bias"] == True:
                 biaslevel = self.bias_level
                 overscan = biaslevel-2
@@ -746,7 +783,7 @@ class Chip(FocalPlane):
                     cr_map[cr_map < 0] = 0
                     crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
                     del cr_map
-                    datetime_obs = datetime.fromtimestamp(timestamp_obs)
+                    datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                     date_obs = datetime_obs.strftime("%y%m%d")
                     time_obs = datetime_obs.strftime("%H%M%S")
                     self.outputCal(
@@ -815,7 +852,7 @@ class Chip(FocalPlane):
                 DarkCombImg.quantize()
                 DarkCombImg = galsim.ImageUS(DarkCombImg)
                 # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
-                datetime_obs = datetime.fromtimestamp(timestamp_obs)
+                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60

ObservationSim/Instrument/Filter.py

@@ -17,8 +17,9 @@ class Filter(object):
         self.filter_id = filter_id
         self.filter_type = filter_type
         self.ccd_bandpass = ccd_bandpass
+        # Load basic filter parameters.
+        # Assume t_exp = 150s for limiting/saturation magnitudes.
         self._getParam(filter_param, filter_type)
-        # if self.filter_dir is not None:
         self.bandpass_full, self.bandpass_sub_list = self._get_bandpasses()
         self.survey_type = self._getSurveyType()
 
@@ -37,7 +38,6 @@ class Filter(object):
         self.sky_background = filter_param.param[filter_type][5]
         self.mag_saturation = filter_param.param[filter_type][6]
         self.mag_limiting = filter_param.param[filter_type][7]
-        # self.filter_dir = filter_param.filter_dir
 
     def is_too_bright(self, mag, margin=-2.5):
         return mag <= self.mag_saturation + margin
@@ -48,9 +48,6 @@ class Filter(object):
 
     def _get_bandpasses(self, filter_dir=None, unit='A'):
         if self.filter_id < 7:  # Photometric
-            # Get full-bandpass
-            # filter_file = os.path.join(filter_dir, self.filter_type+".dat")
-            # bandpass_full = galsim.Bandpass(filter_file, wave_type=unit)
             try:
                 with pkg_resources.files('ObservationSim.Instrument.data.filters').joinpath(self.filter_type.lower() + '.txt') as filter_file:
                     self.filter_bandpass = galsim.Bandpass(str(filter_file), wave_type=unit)
@@ -127,3 +124,5 @@ class Filter(object):
         throughput_file = self.filter_type.lower() + '_throughput.txt'
 
         self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)
+
+        print("for filter %s: mag_limiting: %.3f, mag_saturation: %.3f"%(self.filter_type, self.mag_limiting, self.mag_saturation))

ObservationSim/Instrument/FilterParam.py

@@ -33,6 +33,7 @@ class FilterParam(object):
                     "i":     [7653.2,  1588.1,  6760.0,   8550.0,  0.4960, 0.212, 16.7, 25.9],
                     "z":     [9600.6,  2490.5,  8240.0,  11000.0,  0.2000, 0.123, 15.7, 25.2],
                     "y":     [10051.0, 1590.6,  9130.0,  11000.0,  0.0960, 0.037, 14.4, 24.4],
+                    "FGS":   [5000.0,  8000.0,  3000.0,  11000.0,  0.6500, 0.164, 0., 30.], # [TODO]
                     
                     "GU":    [0.0,        0.0,  2550.0,   4200.0,     1.0, 0.037, 14.0, 26.0],
                     "GV":    [0.0,        0.0,  4000.0,   6500.0,     1.0, 0.037, 14.0, 26.0],

ObservationSim/Instrument/FocalPlane.py

@@ -4,8 +4,8 @@ import numpy as np
 class FocalPlane(object):
     def __init__(self, config=None, survey_type='Photometric', bad_chips=None):
         """Get the focal plane layout
-        Note: (If applicable) make sure the config file have the most basic info, e.g.: "nchip_x"...
         """
+        self.nchips = 42
         if bad_chips == None:
             self.bad_chips = []
         else:
@@ -13,22 +13,24 @@ class FocalPlane(object):
         
         self.ignore_chips = []
         if survey_type == 'Photometric':
-            # for i in range(6):
-            #     self.ignore_chips.append(i+1)
-            #     self.ignore_chips.append(i+25)
             for i in range(5):
                 self.ignore_chips.append(i+1)
                 self.ignore_chips.append(i+26)
             self.ignore_chips.append(10)
             self.ignore_chips.append(21)
+            for i in range(31, 43):
+                self.ignore_chips.append(i)
         elif survey_type == 'Spectroscopic':
-            # for i in range(18):
-            #     self.ignore_chips.append(i+7)
             for i in range(6, 26):
                 if i == 10 or i == 21:
                     continue
                 else:
                     self.ignore_chips.append(i)
+            for i in range(31, 43):
+                self.ignore_chips.append(i)
+        elif survey_type == 'FGS':
+            for i in range(1, 31):
+                self.ignore_chips.append(i)
 
         if config is not None:
             self.nchip_x    = config["nchip_x"]
@@ -42,12 +44,6 @@ class FocalPlane(object):
             if "bad_chips" in config:
                 self.bad_chips = config["bad_chips"]
         else:
-            # self.nchip_x = 5
-            # self.nchip_y = 6
-            # self.npix_tot_x = 49752
-            # self.npix_tot_y = 59516
-            # self.npix_gap_x = 898
-            # self.npix_gap_y = (534, 1309)
             self.nchip_x = 6
             self.nchip_y = 5
             self.npix_tot_x = 59516
@@ -58,8 +54,6 @@ class FocalPlane(object):
         self._getCenter()
 
     def _getCenter(self):
-        # self.cen_pix_x = 0.5 * (self.npix_tot_x + 1.0)
-        # self.cen_pix_y = 0.5 * (self.npix_tot_y + 1.0)
         self.cen_pix_x = 0
         self.cen_pix_y = 0
 

ObservationSim/Instrument/Telescope.py

@@ -42,6 +42,3 @@ class Telescope(object):
 			opticsEff[str(columns[0])] = float(columns[2])
 		f.close()
 		return opticsEff
\ No newline at end of file
-
-	def _get_effCurve(self, effCurve_path):
-		pass
\ No newline at end of file

ObservationSim/Instrument/_util.py

@@ -14,6 +14,17 @@ VC_A = 2.99792458e+18  # speed of light: A/s
 VC_M = 2.99792458e+8   # speed of light: m/s
 H_PLANK = 6.626196e-27 # Plank constant: erg s
 
+def rotate_conterclockwise(x0, y0, x, y, angle):
+	"""
+    Rotate a point counterclockwise by a given angle around a given origin.
+
+    The angle should be given in radians.
+    """
+	angle = np.deg2rad(angle)
+	qx = x0 + np.cos(angle)*(x - x0) - np.sin(angle) * (y - y0)
+	qy = y0 + np.sin(angle)*(x - x0) + np.cos(angle) * (y - y0)
+	return qx, qy
+
 def photonEnergy(lambd):
 	""" The energy of photon at a given wavelength
 
@@ -26,23 +37,22 @@ def photonEnergy(lambd):
 	eph = H_PLANK * nu
 	return eph
 
-'''
-description: 
-param {*} aperture: unit m, default 2 m
-param {*} psf_fwhm: psf fwhm, default 0.1969"
-param {*} pixelSize: pixel size, default 0.074"
-param {*} pmRation: the ratio of souce flux in the limit mag calculation
-param {*} throughputFn: throuput file name
-param {*} readout: unit, e-/pixel
-param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
-param {*} darknoise: unit, e-/pixel/s
-param {*} exTime: exposure time one time, default 150s
-param {*} exNum: exposure number, defautl 1
-param {*} fw, full well value( or saturation value),default 90000e-/pixel
-return {*} limit mag and saturation mag
-'''
-
 def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRation = 0.8, throughputFn = 'i_throughput.txt', readout = 5.0, skyFn= 'sky_emiss_hubble_50_50_A.dat', darknoise = 0.02,exTime = 150, exNum = 1, fw = 90000):
+	'''
+	description: 
+	param {*} aperture: unit m, default 2 m
+	param {*} psf_fwhm: psf fwhm, default 0.1969"
+	param {*} pixelSize: pixel size, default 0.074"
+	param {*} pmRation: the ratio of souce flux in the limit mag calculation
+	param {*} throughputFn: throuput file name
+	param {*} readout: unit, e-/pixel
+	param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
+	param {*} darknoise: unit, e-/pixel/s
+	param {*} exTime: exposure time one time, default 150s
+	param {*} exNum: exposure number, defautl 1
+	param {*} fw, full well value( or saturation value),default 90000e-/pixel
+	return {*} limit mag and saturation mag
+	'''
 	try:
 		with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(throughputFn) as data_file:
 			throughput_f = np.loadtxt(data_file)