C6 debugging

Catalog/C3Catalog.py

@@ -105,7 +105,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(
@@ -138,8 +138,8 @@ class C3Catalog(CatalogBase):
                 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
@@ -199,7 +199,7 @@ class C3Catalog(CatalogBase):
             pmdec_list = np.zeros(nstars).tolist()
             rv_list = np.zeros(nstars).tolist()
             parallax_list = [1e-9] * nstars
-            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(

Catalog/C6_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
+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
+
+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 C6_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.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
+
+        # (DEBUG)
+        self.max_size = 0.
+
+        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"] and not config["run_option"]["galaxy_only"]:
+            star_file = config["input_path"]["star_cat"]
+            star_SED_file = config["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["input_path"] and config["input_path"]["galaxy_cat"] and not config["run_option"]["star_only"]:
+            galaxy_dir = config["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["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+        
+        if "rotateEll" in config["shear_setting"]:
+            self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
+        else:
+            self.rotation = 0.
+
+        # 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)
+        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_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_true_g_lsst'][igals]
+            # 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['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'] = str(pix_id) + '%02d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+            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
+        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 > 10:
+            #     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'] = self.ids
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # 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):
+        self.objs = []
+        self.ids = 0
+        
+        if "star_cat" in self.config["input_path"] and self.config["input_path"]["star_cat"] and not self.config["run_option"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['catalog']
+            for pix in self.pix_list:
+                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["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
+            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:
+                        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':
+            dist_L_pc = (1 + obj.z) * comoving_dist(z=obj.z)[0]
+            factor = (10 / dist_L_pc)**2
+            flux = np.matmul(self.pcs, obj.coeff) * factor
+            if np.any(flux < 0):
+                raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+            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]
+            # 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 == 'galaxy' or obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            # print("sed_photon = ", sed_photon)
+            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)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+        
+        del wave
+        del flux
+        return sed

Catalog/NGPCatalog.py

@@ -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(
@@ -144,8 +144,8 @@ class NGPCatalog(CatalogBase):
             #     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(

ObservationSim/Config/ChipOutput.py

@@ -40,18 +40,25 @@ class ChipOutput(object):
         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"
+        # 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"
 
-        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
 
-        self.hdr = hdr1 + hdr2 + hdr3 + hdr4
-        self.fmt = fmt1 + fmt2 + fmt3 + fmt4
+        hdr1  = "# obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type "
+        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))
 
@@ -78,8 +85,8 @@ class ChipOutput(object):
         
         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.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/Instrument/Chip/Chip.py

@@ -419,7 +419,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(
@@ -585,7 +585,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
@@ -689,7 +689,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
@@ -746,7 +746,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 +815,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/MockObject/CatalogBase.py

@@ -38,8 +38,8 @@ class CatalogBase(metaclass=ABCMeta):
             "mag_use_normal":100.,
             "theta":0.,
             "kappa":0.,
-            "gamma1":0.,
-            "gamma2":0.,
+            "g1":0.,
+            "g2":0.,
             "bfrac":0.,
             "av":0.,
             "redden":0.,
@@ -48,11 +48,24 @@ class CatalogBase(metaclass=ABCMeta):
             "ell_bulge":0.,
             "ell_disk":0.,
             "ell_tot":0.,
+            "e1_disk":0.,
+            "e2_disk":0.,
+            "e1_bulge":0.,
+            "e2_bulge":0.,
             "teff":0.,
             "logg":0.,
             "feh":0.,
-            "g1":0.,
-            "g2":0.,
+            # C6 galaxies parameters
+            "e1":0.,
+            "e2":0.,
+            "bulgemass":0.,
+            "diskmass":0.,
+            "size":0.,
+            "detA":0.,
+            "type":0,
+            "veldisp":0.,
+            "coeff": np.zeros(20),
+            # Astrometry related
             "pmra":0.,
             "pmdec":0.,
             "rv":0.,

ObservationSim/MockObject/Galaxy.py

@@ -12,18 +12,18 @@ from ObservationSim.MockObject.MockObject import MockObject
 class Galaxy(MockObject):
     def __init__(self, param, rotation=None, logger=None):
         super().__init__(param, logger=logger)
-        self.thetaR = self.param["theta"]
-        self.bfrac = self.param["bfrac"]
-        self.hlr_disk = self.param["hlr_disk"]
-        self.hlr_bulge = self.param["hlr_bulge"]
+        # self.thetaR = self.param["theta"]
+        # self.bfrac = self.param["bfrac"]
+        # self.hlr_disk = self.param["hlr_disk"]
+        # self.hlr_bulge = self.param["hlr_bulge"]
 
         # Extract ellipticity components
-        self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees)
-        self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees)
-        self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees)
-        self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2
-        self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2
-        self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2
+        # self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees)
+        # self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees)
+        # self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees)
+        # self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2
+        # self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2
+        # self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2
 
         if rotation is not None:
             self.rotateEllipticity(rotation)
@@ -147,12 +147,12 @@ class Galaxy(MockObject):
             bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
             bulge = bulge.shear(bulge_shape)
 
-            # gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
+            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
 
             # (TEST) Random knots
-            knots = galsim.RandomKnots(npoints=100, profile=disk)
-            kfrac = np.random.random()*(1.0 - self.bfrac)
-            gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
+            # knots = galsim.RandomKnots(npoints=100, profile=disk)
+            # kfrac = np.random.random()*(1.0 - self.bfrac)
+            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
 
             gal = gal.withFlux(nphotons)
             gal_shear = galsim.Shear(g1=g1, g2=g2)
@@ -287,12 +287,12 @@ class Galaxy(MockObject):
             bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
             bulge = bulge.shear(bulge_shape)
 
-            # gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
+            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
 
             # (TEST) Random knots
-            knots = galsim.RandomKnots(npoints=100, profile=disk)
-            kfrac = np.random.random()*(1.0 - self.bfrac)
-            gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
+            # knots = galsim.RandomKnots(npoints=100, profile=disk)
+            # kfrac = np.random.random()*(1.0 - self.bfrac)
+            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
 
             gal = gal.withFlux(tel.pupil_area * exptime)
             gal_shear = galsim.Shear(g1=g1, g2=g2)

ObservationSim/MockObject/MockObject.py

@@ -12,7 +12,10 @@ from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 class MockObject(object):
     def __init__(self, param, logger=None):
         self.param = param
-
+        
+        for key in self.param:
+            setattr(self, key, self.param[key])
+        
         if self.param["star"] == 0:
             self.type = "galaxy"
         elif self.param["star"] == 1:
@@ -20,32 +23,32 @@ class MockObject(object):
         elif self.param["star"] == 2:
             self.type = "quasar"
 
-        self.id = self.param["id"]
-        self.ra = self.param["ra"]
-        self.dec = self.param["dec"]
-        self.ra_orig = self.param["ra_orig"]
-        self.dec_orig = self.param["dec_orig"]
-        self.z = self.param["z"]
-        self.sed_type = self.param["sed_type"]
-        self.model_tag = self.param["model_tag"]
-        self.mag_use_normal = self.param["mag_use_normal"]
+        # self.id = self.param["id"]
+        # self.ra = self.param["ra"]
+        # self.dec = self.param["dec"]
+        # self.ra_orig = self.param["ra_orig"]
+        # self.dec_orig = self.param["dec_orig"]
+        # self.z = self.param["z"]
+        # self.sed_type = self.param["sed_type"]
+        # self.model_tag = self.param["model_tag"]
+        # self.mag_use_normal = self.param["mag_use_normal"]
         self.sed = None
 
         # Place holder for outputs
-        self.av = self.param["av"]
-        self.redden = self.param["redden"]
-        self.pmra = self.param["pmra"]
-        self.pmdec = self.param["pmdec"]
-        self.rv = self.param["rv"]
-        self.parallax = self.param["parallax"]
-        self.g1 = self.param["g1"]
-        self.g2 = self.param["g2"]
-        self.thetaR = self.param["theta"]
-        self.bfrac = self.param["bfrac"]
-        self.hlr_disk = self.param["hlr_disk"]
-        self.hlr_bulge = self.param["hlr_bulge"]
-        self.e1_disk, self.e2_disk = 0., 0.
-        self.e1_bulge, self.e2_bulge = 0., 0.
+        # self.av = self.param["av"]
+        # self.redden = self.param["redden"]
+        # self.pmra = self.param["pmra"]
+        # self.pmdec = self.param["pmdec"]
+        # self.rv = self.param["rv"]
+        # self.parallax = self.param["parallax"]
+        # self.g1 = self.param["g1"]
+        # self.g2 = self.param["g2"]
+        # self.thetaR = self.param["theta"]
+        # self.bfrac = self.param["bfrac"]
+        # self.hlr_disk = self.param["hlr_disk"]
+        # self.hlr_bulge = self.param["hlr_bulge"]
+        # self.e1_disk, self.e2_disk = 0., 0.
+        # self.e1_bulge, self.e2_bulge = 0., 0.
         self.additional_output_str = ""
 
         self.logger = logger

ObservationSim/MockObject/_util.py

 import numpy as np
 import os
 from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
-from scipy import interpolate
+from scipy import interpolate, integrate
 from astropy.table import Table
 import galsim
 
@@ -9,6 +9,14 @@ 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 comoving_dist(z, om_m=0.3111, om_L=0.6889, h=0.6766): 
+    # Return comving distance in pc
+    H0 = h*100. # km / (s Mpc)
+    def dist_int(z):
+        return 1./np.sqrt(om_m*(1.+z)**3 + om_L)
+    res, err = integrate.quad(dist_int, 0., z)
+    return [res * (VC_M/1e3/H0) * 1e6, err * (VC_M/1e3/H0) * 1e6]
+
 def magToFlux(mag):
     """
     flux of a given AB magnitude

ObservationSim/ObservationSim.py

@@ -61,7 +61,7 @@ class Observation(object):
 
         # print(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
         # print("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
-        # print("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat())
+        # print("Time: %s" % datetime.utcfromtimestamp(pointing.timestamp).isoformat())
         # print("Exposure time: %f" % pointing.exp_time)
         # print("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
         # print("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
@@ -70,7 +70,7 @@ class Observation(object):
         # print(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
         chip_output.logger.info(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
         chip_output.logger.info("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
-        chip_output.logger.info("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat())
+        chip_output.logger.info("Time: %s" % datetime.utcfromtimestamp(pointing.timestamp).isoformat())
         chip_output.logger.info("Exposure time: %f" % pointing.exp_time)
         chip_output.logger.info("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
         chip_output.logger.info("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
@@ -92,7 +92,7 @@ class Observation(object):
 
         # Apply astrometric simulation for pointing
         if self.config["obs_setting"]["enable_astrometric_model"]:
-            dt = datetime.fromtimestamp(pointing.timestamp)
+            dt = datetime.utcfromtimestamp(pointing.timestamp)
             date_str = dt.date().isoformat()
             time_str = dt.time().isoformat()
             ra_cen, dec_cen = on_orbit_obs_position(
@@ -163,7 +163,7 @@ class Observation(object):
 
         if pointing.pointing_type == 'MS':
             # Load catalogues and templates
-            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output)
+            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt)
             chip_output.create_output_file()
             self.nobj = len(self.cat.objs)
 
@@ -226,10 +226,13 @@ class Observation(object):
                     mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()],
                     margin=self.config["obs_setting"]["mag_sat_margin"]):
                     # print("obj too birght!!", flush=True)
-                    if obj.type != 'galaxy':
-                        bright_obj += 1
-                        obj.unload_SED()
-                        continue
+                    # if obj.type != 'galaxy':
+                    #     bright_obj += 1
+                    #     obj.unload_SED()
+                    #     continue
+                    bright_obj += 1
+                    obj.unload_SED()
+                    continue
                 if filt.is_too_dim(
                     mag=obj.getMagFilter(filt),
                     margin=self.config["obs_setting"]["mag_lim_margin"]):
@@ -366,7 +369,7 @@ class Observation(object):
                 logger=chip_output.logger)
             
             if pointing.pointing_type == 'MS':
-                datetime_obs = datetime.fromtimestamp(pointing.timestamp)
+                datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
                 date_obs = datetime_obs.strftime("%y%m%d")
                 time_obs = datetime_obs.strftime("%H%M%S")
                 h_prim = generatePrimaryHeader(

ObservationSim/_util.py

@@ -144,7 +144,7 @@ def makeSubDir_PointingList(path_dict, config, pointing_ID=0):
     return subImgdir, prefix
 
 def get_shear_field(config, shear_cat_file=None):
-    if not config["shear_setting"]["shear_type"] in ["constant", "extra"]:
+    if not config["shear_setting"]["shear_type"] in ["constant", "extra", "catalog"]:
         raise ValueError("Please set a right 'shear_method' parameter.")
 
     if config["shear_setting"]["shear_type"] == "constant":
@@ -152,16 +152,19 @@ def get_shear_field(config, shear_cat_file=None):
         g2 = config["shear_setting"]["reduced_g2"]
         nshear = 1
         # TODO logging
-    else:
+    elif config["shear_setting"]["shear_type"] == "extra":
         # TODO logging
         if not os.path.exists(shear_cat_file):
-            raise ValueError("Cannot find shear catalog file.")
+            raise ValueError("Cannot find external shear catalog file.")
         try:
             shearCat = np.loadtxt(shear_cat_file)
             nshear = shearCat.shape[0]
             g1, g2 = shearCat[:, 0], shearCat[:, 1]
         except:
-            print("Failed to the shear catalog file.")
+            print("Failed to read the shear catalog file.")
             print("Setting to no shear.")
             g1, g2 = 0., 0.
+    else:
+        g1, g2 = 0., 0.
+        nshear = 0
     return g1, g2, nshear
\ No newline at end of file

config/config_NGP_dev.yaml

@@ -9,17 +9,17 @@
 # Base diretories and naming setup
 # Can add some of the command-line arguments here as well;
 # OK to pass either way or both, as long as they are consistent
-# work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/"
-work_dir: "/share/"
+work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/"
+# work_dir: "/share/"
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "10-pointings_for_Tianmeng"
+run_name: "test_mag_completeness_no_extinction-reddening"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
 # - pointing_dir: null
 # - pointing_file: null
-# pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-pointing_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/"
+pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/"
+# pointing_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/"
 pointing_file: "pointing_test_NGP_2.17.dat"
 
 # Whether to use MPI
@@ -32,7 +32,7 @@ run_option:
 
   # Output catalog only?
   # If yes, no imaging simulation will run
-  out_cat_only: NO
+  out_cat_only: YES
 
   # Only simulate stars?
   star_only: NO
@@ -74,7 +74,7 @@ obs_setting:
   # - give a list of indexes of pointings: [ip_1, ip_2...]
   # - run all pointings: null
   # Note: only valid when a pointing list is specified
-  run_pointings: null
+  run_pointings: [0, 1, 2]
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
@@ -87,7 +87,7 @@ obs_setting:
   enable_astrometric_model: True
 
   # Cut by saturation/limiting magnitude in which band?
-  cut_in_band: "g"
+  cut_in_band: "z"
 
   # saturation magnitude margin
   mag_sat_margin: -2.5
@@ -102,10 +102,10 @@ obs_setting:
 # Default path settings for WIDE survey simulation
 input_path:
   cat_dir: "OnOrbitCalibration/CTargets20211231"
-  # star_cat: "CT-NGP_r1.8_G28.hdf5"
-  # galaxy_cat: "galaxyCats_r_3.0_healpix_shift_192.859500_27.128300.hdf5"
-  star_cat: "C4_MMW_Cluster_D50_GGC_Astrometry_G.hdf5"
-  galaxy_cat: "galaxyCats_r_10.0_healpix.hdf5"
+  star_cat: "CT-NGP_r1.8_G28.hdf5"
+  galaxy_cat: "galaxyCats_r_3.0_healpix_shift_192.859500_27.128300.hdf5"
+  # star_cat: "C4_MMW_Cluster_D50_GGC_Astrometry_G.hdf5"
+  # galaxy_cat: "galaxyCats_r_10.0_healpix.hdf5"
 
 SED_templates_path:
   star_SED: "Catalog_20210126/SpecLib.hdf5"

run_C6.pbs

+#!/bin/bash
+
+#PBS -N SIMS
+#PBS -l nodes=wcl-1:ppn=16+wcl-2:ppn=16+wcl-3:ppn=16+wcl-4:ppn=16+wcl-5:ppn=16+wcl-6:ppn=16
+#PBS -u fangyuedong
+###PBS -j oe
+
+cd $PBS_O_WORKDIR
+
+NP=96
+date
+
+mpirun -np $NP python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_C6.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config
+

run_sim.py

@@ -100,7 +100,9 @@ if __name__=='__main__':
     # run_sim(Catalog=C3Catalog)
 
     # To run calibration field NGP simulation
-    from Catalog.NGPCatalog import NGPCatalog
-    run_sim(Catalog=NGPCatalog)
+    # from Catalog.NGPCatalog import NGPCatalog
+    # run_sim(Catalog=NGPCatalog)
     # from Catalog.NJU_Catalog import NJU_Catalog
     # run_sim(Catalog=NJU_Catalog)
+    from Catalog.C6_Catalog import C6_Catalog
+    run_sim(Catalog=C6_Catalog)

test_C6.sh

+#!/bin/bash
+
+date
+
+# mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/test/CSST/config
+# python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config
+python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_C6.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config
+