update output catalog, apply astrometry module to galaxy catalog, add new PSF tests

Catalog/C3Catalog.py

@@ -84,10 +84,46 @@ class C3Catalog(CatalogBase):
         self.rng_sedGal = random.Random()
         self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
         self.ud = galsim.UniformDeviate(pix_id)
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra_true'][:]
+        dec_arr = gals['dec_true'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.fromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=ngals,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2015.5",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
         for igals in range(ngals):
             param = self.initialize_param()
-            param['ra'] = gals['ra_true'][igals]
-            param['dec'] = gals['dec_true'][igals]
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra_true'][igals]
+            param['dec_orig'] = gals['dec_true'][igals]
             if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
                 continue
             param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
@@ -129,7 +165,8 @@ class C3Catalog(CatalogBase):
                 param['star'] = 2 # Quasar
 
             self.ids += 1
-            param['id'] = self.ids
+            # param['id'] = self.ids
+            param['id'] = gals['galaxyID'][igals]
             
             if param['star'] == 0:
                 obj = Galaxy(param, self.rotation)
@@ -141,9 +178,9 @@ class C3Catalog(CatalogBase):
     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"][:]
-        # if "astrometric_lib" in self.config["obs_setting"] and self.config["obs_setting"]["enable_astrometric_model"]:
         if self.config["obs_setting"]["enable_astrometric_model"]:
             ra_list = ra_arr.tolist()
             dec_list = dec_arr.tolist()
@@ -170,20 +207,22 @@ class C3Catalog(CatalogBase):
                 input_vz=self.pointing.sat_vz,
                 input_epoch="J2015.5",
                 input_date_str=date_str,
-                input_time_str=time_str,
-                # lib_path=lib_path
+                input_time_str=time_str
             )
         for istars in range(nstars):
             param = self.initialize_param()
             param['ra'] = ra_arr[istars]
             param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
             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'] = 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]

ObservationSim/Config/ChipOutput.py

@@ -23,73 +23,81 @@ class ChipOutput(object):
         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.subdir = subdir
 
-        # hdr1  = "#ID ID_chip filter xImage yImage ra dec z mag flag SNR "
-        hdr1  = "#ID ID_chip filter xImage yImage ra dec z mag flag "
-        hdr2  = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge "
-        hdr3  = "e1PSF e2PSF e1 e2 g1 g2 e1OBS e2OBS "
-        hdr4  = "sed_type av redden "
-        hdr5  = "star_model teff logg feh\n"
-        # fmt1  = "%10d %4d %5s %10.3f %10.3f %15.6f %15.6f %7.4f %8.4f %2d %9.2f "
-        fmt1  = "%10d %4d %5s %10.3f %10.3f %15.6f %15.6f %7.4f %8.4f %2d "
-        fmt2  = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f "
-        fmt3  = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f "
-        fmt4 = "%2d %8.4f %8.4f "
-        fmt5 = "%10s %8.4f %8.4f %8.4f\n"
-        self.hdr = hdr1 + hdr2 + hdr3 + hdr4 + hdr5
-        self.fmt = fmt1 + fmt2 + fmt3 + fmt4 + fmt5
+        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\n"
+
+        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\n"
+
+        self.hdr = hdr1 + hdr2 + hdr3 + hdr4
+        self.fmt = fmt1 + fmt2 + fmt3 + fmt4
+
         print("pointing_type = %s\n"%(pointing_type))
         if pointing_type == 'MS':
             self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
             print("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
             self.cat.write(self.hdr)
 
-    def updateHDR(self, hdr):
-        hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise,        "comment":"read noise in e-/pixel"},
-                {"name":"DARK",    "value":self.chip.dark_noise,        "comment":"Dark noise (e-/pixel/s)"},
-                {"name":"EXPTIME", "value":self.exptime,          "comment":"exposure time in second"},
-                {"name":"GAIN",    "value":self.chip.gain,             "comment":"CCD gain in e-/ADU"},
-                {"name":"SATURATE","value":65535.0,          "comment":"saturation level"},
-                {"name":"CCDCHIP",  "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"},
-                {"name":"FILTER",  "value":self.filt.filter_type,          "comment":"filter name"},
-                {"name":"MJD-OBS", "value":self.mjdTime,          "comment":"Modified Julian Date (MJD) of observation"},
-                {"name":"DATE-OBS","value":self.imgKey1,          "comment":"Date of observation"},
-                {"name":"EQUINOX", "value":2000.0},
-                {"name":"RADECSYS","value":"ICRS"},
-                {"name":"RA",      "value":self.ra_cen,           "comment":"telescope pointing center"},
-                {"name":"DEC",     "value":self.dec_cen,          "comment":"telescope pointing center"},
-                {"name":"OBJECT",  "value":"CSS-OS"},
-                {"name":"WCSDIM",  "value":2.0,              "comment":"WCS Dimensionality"},
-                {"name":"EXTNAME", "value":"IM1",            "comment":"Extension name"},
-                {"name":"BSCALE",  "value":1.0},
-                {"name":"BZERO",   "value":0.0},
-                {"name":"OBSID",   "value":self.imgKey0,          "comment":"Observation ID"},
-                {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"},
-                {"name":"RSPEED",  "value":10.0,             "comment":"Read speed"},
-                {"name":"CHIPTEMP","value":-100.0,           "comment":"Chip temperature"},
-                {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"},
-                {"name":"CCDSUM",  "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
-                {"name":"NSUM",    "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
-                {"name":"AUTHOR",  "value":"CSST-Sim Group"},
-                {"name":"GROUP",   "value":"Weak Lensing Working Group for CSST"}]
-        for item in hdrNew:
-            hdr.add_record(item)
-        return hdr
+    # def updateHDR(self, hdr):
+    #     hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise,        "comment":"read noise in e-/pixel"},
+    #             {"name":"DARK",    "value":self.chip.dark_noise,        "comment":"Dark noise (e-/pixel/s)"},
+    #             {"name":"EXPTIME", "value":self.exptime,          "comment":"exposure time in second"},
+    #             {"name":"GAIN",    "value":self.chip.gain,             "comment":"CCD gain in e-/ADU"},
+    #             {"name":"SATURATE","value":65535.0,          "comment":"saturation level"},
+    #             {"name":"CCDCHIP",  "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"},
+    #             {"name":"FILTER",  "value":self.filt.filter_type,          "comment":"filter name"},
+    #             {"name":"MJD-OBS", "value":self.mjdTime,          "comment":"Modified Julian Date (MJD) of observation"},
+    #             {"name":"DATE-OBS","value":self.imgKey1,          "comment":"Date of observation"},
+    #             {"name":"EQUINOX", "value":2000.0},
+    #             {"name":"RADECSYS","value":"ICRS"},
+    #             {"name":"RA",      "value":self.ra_cen,           "comment":"telescope pointing center"},
+    #             {"name":"DEC",     "value":self.dec_cen,          "comment":"telescope pointing center"},
+    #             {"name":"OBJECT",  "value":"CSS-OS"},
+    #             {"name":"WCSDIM",  "value":2.0,              "comment":"WCS Dimensionality"},
+    #             {"name":"EXTNAME", "value":"IM1",            "comment":"Extension name"},
+    #             {"name":"BSCALE",  "value":1.0},
+    #             {"name":"BZERO",   "value":0.0},
+    #             {"name":"OBSID",   "value":self.imgKey0,          "comment":"Observation ID"},
+    #             {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"},
+    #             {"name":"RSPEED",  "value":10.0,             "comment":"Read speed"},
+    #             {"name":"CHIPTEMP","value":-100.0,           "comment":"Chip temperature"},
+    #             {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"},
+    #             {"name":"CCDSUM",  "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
+    #             {"name":"NSUM",    "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
+    #             {"name":"AUTHOR",  "value":"CSST-Sim Group"},
+    #             {"name":"GROUP",   "value":"Weak Lensing Working Group for CSST"}]
+    #     for item in hdrNew:
+    #         hdr.add_record(item)
+    #     return hdr
 
     # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2):
-    def cat_add_obj(self, obj, pos_img, pos_shear, g1, g2):
+    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
-        e1, e2, g1, g2, e1OBS, e2OBS = obj.getObservedEll(g1, g2)
-        if obj.type == 'galaxy':
-            line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge,
-                obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge,
-                pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0)
-        elif obj.type == "quasar":
-            line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z,
-                               obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-                               pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0)
-        else:
-            line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
-                pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh'])
-        # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS)
+        # if obj.type == 'galaxy':
+        #     line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge,
+        #         obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge,
+        #         pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0)
+        # elif obj.type == "quasar":
+        #     line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z,
+        #                        obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        #                        pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0)
+        # else:
+        #     line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+        #         pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh'])
+        # print(
+        #     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 = 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)
         self.cat.write(line)

ObservationSim/MockObject/Catalog.py

-import os
-import numpy as np
-import random
-import galsim
-import h5py as h5
-import healpy as hp
-from astropy.table import Table
-from astropy.coordinates import spherical_to_cartesian
-
-from ObservationSim.MockObject.Star import Star
-from ObservationSim.MockObject.Galaxy import Galaxy
-from ObservationSim.MockObject.Quasar import Quasar
-from ObservationSim.MockObject._util import seds, sed_assign, extAv
-
-NSIDE = 128
-
-class Catalog(object):
-    def __init__(self, config, chip, cat_dir=None, pRa=None, pDec=None, sed_dir=None, rotation=None):
-        if cat_dir is not None:
-            self.cat_dir = cat_dir
-        else:
-            self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
-        self.sed_dir = sed_dir
-        self.chip = chip
-        
-        self.seed_Av = config["random_seeds"]["seed_Av"]
-        
-        if pRa is not None:
-            self.pRa = float('{:8.4f}'.format(pRa))
-        if pDec is not None:
-            self.pDec = float('{:8.4f}'.format(pDec))
-
-        # star_file = 'stars_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5'
-        # galaxy_file = 'galaxies_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5'
-
-        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"]
-            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"]:
-            galaxy_file = config["input_path"]["galaxy_cat"]
-            self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
-            self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
-            self._load_SED_lib_gals()
-        if "rotateEll" in config["shear_setting"]:
-            self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
-        else:
-            self.rotation = 0.
-        # self._load_SED_info()
-        self._get_healpix_list()
-        self._load()
-
-    def _get_healpix_list(self):
-        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
-        # self.sky_coverage_enlarged = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.5)
-        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]))
-        # phi, theta = ra, np.pi/2. - dec
-        # vertices_pix = np.unique(hp.ang2pix(NSIDE, theta, phi))
-        # print(vertices_pix)
-        vertices = spherical_to_cartesian(1., dec, ra)
-        self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
-        print("HEALPix List: ", self.pix_list)
-
-    def _load_SED_lib_star(self):
-        self.tempSED_star = h5.File(self.star_SED_path,'r')
-
-    def _load_SED_lib_gals(self):
-        # tempdir_gal = os.path.join(self.template_dir, "Galaxy/")
-        # tempdir_star = os.path.join(self.template_dir, "PicklesStars/")
-        # self.tempSed_star, self.tempRed_star = seds("star.list", seddir=tempdir_star)
-        self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path)
-
-    def _load_SED_info(self):
-        sed_info_file = os.path.join(self.sed_dir, "sed.info")
-        sed_info = Table.read(sed_info_file, format="ascii")
-        self.cosids = sed_info["IDCosmoDC2"]
-        self.objtypes = sed_info["objType"]
-
-    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)
-        for igals in range(ngals):
-            param = {}
-            param['ra'] = gals['ra_true'][igals]
-            param['dec'] = gals['dec_true'][igals]
-            param['z'] = gals['redshift_true'][igals]
-            param['model_tag'] = 'None'
-            param['gamma1'] = 0
-            param['gamma2'] = 0
-            param['kappa'] = 0
-            param['delta_ra'] = 0
-            param['delta_dec'] = 0
-            sersicB = gals['sersic_bulge'][igals]
-            hlrMajB = gals['size_bulge_true'][igals]
-            hlrMinB = gals['size_minor_bulge_true'][igals]
-            sersicD = gals['sersic_disk'][igals]
-            hlrMajD = gals['size_disk_true'][igals]
-            hlrMinD = gals['size_minor_disk_true'][igals]
-            aGal = gals['size_true'][igals]
-            bGal = gals['size_minor_true'][igals]
-            param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
-            param['theta'] = gals['position_angle_true'][igals]
-            param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
-            param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
-            param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
-            param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
-            param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
-
-            # Assign each galaxy a template SED
-            param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
-            param['redden'] = self.tempRed_gal[param['sed_type']]
-            param['av'] = self.avGal[int(self.ud()*self.nav)]
-            if param['sed_type'] <= 5:
-                param['av'] = 0.0
-                param['redden'] = 0
-            param['star'] = 0   # Galaxy
-            if param['sed_type'] >= 29:
-                param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
-                param['star'] = 2 # Quasar
-
-            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
-                continue
-            param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
-            if param['mag_use_normal'] >= 26.5:
-                continue
-            self.ids += 1
-            param['id'] = self.ids
-            
-            if param['star'] == 0:
-                obj = Galaxy(param, self.rotation)
-                self.objs.append(obj)
-            if param['star'] == 2:
-                obj = Quasar(param)
-                self.objs.append(obj)
-
-    def _load_stars(self, stars, pix_id=None):
-        nstars = len(stars['sourceID'])
-        for istars in range(nstars):
-            param = {}
-            param['ra'] = stars['RA'][istars]
-            param['dec'] = stars['Dec'][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['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)
-            self.objs.append(obj)
-
-    def _load(self):
-        self.nav = 15005
-        self.avGal = extAv(self.nav, seed=self.seed_Av)
-        gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
-        star_cat = h5.File(self.star_path, 'r')['catalog']
-        self.objs = []
-        self.ids = 0
-        for pix in self.pix_list:
-            gals = gals_cat[str(pix)]
-            stars = star_cat[str(pix)]
-            self._load_gals(gals, pix_id=pix)
-            self._load_stars(stars, pix_id=pix)
-        print("number of objects in catalog: ", len(self.objs))
-        del self.avGal
\ No newline at end of file

ObservationSim/MockObject/CatalogBase.py

@@ -29,29 +29,33 @@ class CatalogBase(metaclass=ABCMeta):
             "star":-1,
             "id":-1,
             "ra":0,
-            "dec":0,
-            "z":0,
+            "dec":0.,
+            "ra_orig":0.,
+            "dec_orig":0.,
+            "z":0.,
             "sed_type":-1,
             "model_tag":"unknown",
-            "mag_use_normal":100,
-            "theta":0,
-            "kappa":0,
-            "gamma1":0,
-            "gamma2":0,
-            "bfrac":0,
-            "hlr_bulge":0,
-            "hlr_disk":0,
-            "ell_bulge":0,
-            "ell_disk":0,
-            "ell_tot":0,
-            "teff":0,
-            "logg":0,
-            "feh":0,
-            "g1":0,
-            "g2":0,
-            "pmra":0,
-            "pmdec":0,
-            "rv":0,
+            "mag_use_normal":100.,
+            "theta":0.,
+            "kappa":0.,
+            "gamma1":0.,
+            "gamma2":0.,
+            "bfrac":0.,
+            "av":0.,
+            "redden":0.,
+            "hlr_bulge":0.,
+            "hlr_disk":0.,
+            "ell_bulge":0.,
+            "ell_disk":0.,
+            "ell_tot":0.,
+            "teff":0.,
+            "logg":0.,
+            "feh":0.,
+            "g1":0.,
+            "g2":0.,
+            "pmra":0.,
+            "pmdec":0.,
+            "rv":0.,
             "parallax":1e-9
         }
         return param

ObservationSim/MockObject/Galaxy.py

@@ -131,12 +131,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
+            # # (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
 
             gal = gal.withFlux(nphotons)
             gal_shear = galsim.Shear(g1=g1, g2=g2)
@@ -220,12 +220,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

@@ -17,15 +17,34 @@ class MockObject(object):
             self.type = "star"
         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.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.
+
     def getMagFilter(self, filt):
         if filt.filter_type in ["GI", "GV", "GU"]:
             return self.param["mag_use_normal"]
@@ -304,5 +323,5 @@ class MockObject(object):
         snr_obj = img_flux / sig_obj
         return snr_obj
 
-    def getObservedEll(self, g1=0, g2=0):
-        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+    # def getObservedEll(self, g1=0, g2=0):
+    #     return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

ObservationSim/MockObject/__init__.py

@@ -3,6 +3,5 @@ from .Galaxy import Galaxy
 from .CatalogBase import CatalogBase
 from .Quasar import Quasar
 from .Star import Star
-from .Catalog import Catalog
 from .SkybackgroundMap import *
 # from .CosmicRay import CosmicRay
\ No newline at end of file

ObservationSim/ObservationSim.py

@@ -19,10 +19,8 @@ class Observation(object):
     def __init__(self, config, Catalog, work_dir=None, data_dir=None):
         self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir)
         self.config = config
-        # self.tel = Telescope(optEffCurve_path=self.path_dict["mirror_file"]) 
         self.tel = Telescope()
         self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) 
-        # self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) 
         self.filter_param = FilterParam() 
         self.chip_list = []
         self.filter_list = []
@@ -44,10 +42,6 @@ class Observation(object):
             # Make Chip & Filter lists
             chip = Chip(
                 chipID=chipID, 
-                # ccdEffCurve_dir=self.path_dict["ccd_dir"], 
-                # CRdata_dir=self.path_dict["CRdata_dir"], 
-                # normalize_dir=self.path_dict["normalize_dir"], 
-                # sls_dir=self.path_dict["sls_dir"], 
                 config=self.config)
             filter_id, filter_type = chip.getChipFilter()
             filt = Filter(filter_id=filter_id, 
@@ -114,7 +108,6 @@ class Observation(object):
                 skyMap=flat_normal, 
                 blueLimit=filt.blue_limit, 
                 redLimit=filt.red_limit, 
-                # skyfn=self.path_dict["sky_file"], 
                 conf=chip.sls_conf, 
                 pixelSize=chip.pix_scale, 
                 isAlongY=0)
@@ -130,7 +123,7 @@ class Observation(object):
             bright_obj = 0
             dim_obj = 0
             for j in range(self.nobj):
-                # if j >= 100:
+                # if j >= 10:
                 #     break
                 obj = self.cat.objs[j]
                 if obj.type == 'star' and self.config["run_option"]["galaxy_only"]:
@@ -170,13 +163,13 @@ class Observation(object):
 
                 if self.config["shear_setting"]["shear_type"] == "constant":
                     if obj.type == 'star':
-                        obj.param["g1"], obj.param["g2"] = 0, 0
+                        obj.g1, obj.g2 = 0., 0.
                     else:
-                        obj.param["g1"], obj.param["g2"] = self.g1_field, self.g2_field
+                        obj.g1, obj.g2 = self.g1_field, self.g2_field
                 elif self.config["shear_setting"]["shear_type"] == "extra":
                     try:
                         # TODO: every object with individual shear from input catalog(s)
-                        obj.param["g1"], obj.param["g2"] = self.g1_field[j], self.g2_field[j]
+                        obj.g1, obj.g2 = self.g1_field[j], self.g2_field[j]
                     except:
                         print("failed to load external shear.")
                         pass
@@ -205,8 +198,8 @@ class Observation(object):
                             bandpass_list=filt.bandpass_sub_list, 
                             filt=filt, 
                             chip=chip, 
-                            g1=obj.param["g1"], 
-                            g2=obj.param["g2"], 
+                            g1=obj.g1, 
+                            g2=obj.g2, 
                             exptime=pointing.exp_time
                             )
                     elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]:
@@ -217,14 +210,15 @@ class Observation(object):
                             bandpass_list=filt.bandpass_sub_list, 
                             filt=filt, 
                             chip=chip, 
-                            g1=obj.param["g1"], 
-                            g2=obj.param["g2"], 
+                            g1=obj.g1, 
+                            g2=obj.g2, 
                             exptime=pointing.exp_time,
                             normFilter=norm_filt,
                             )
                     if isUpdated:
                         # TODO: add up stats
-                        chip_output.cat_add_obj(obj, pos_img, pos_shear, obj.param["g1"], obj.param["g2"])
+                        # print("updating output catalog...")
+                        chip_output.cat_add_obj(obj, pos_img, pos_shear)
                         pass
                     else:
                         # print("object omitted", flush=True)

config/config_C3.yaml

@@ -67,7 +67,6 @@ obs_setting:
   image_rot: -113.4333
 
   # Number of calibration pointings
-  # Note: only valid when a pointing list is specified
   np_cal: 0
 
   # Run specific pointing(s):
@@ -100,7 +99,6 @@ SED_templates_path:
   star_SED: "Catalog_20210126/SpecLib.hdf5"
   galaxy_SED: "Templates/Galaxy/"
 
-
 ###############################################
 # PSF setting
 ###############################################

config/config_NGP.yaml

+---
+###############################################
+#
+#  Configuration file for CSST simulation
+#  CSST-Sim Group, 2021/10/07, version 0.3
+#
+###############################################
+
+# 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: "/public/home/fangyuedong/sim_code_release/CSST/test/"
+work_dir: "/public/home/fangyuedong/temp/CSST/workplace/"
+data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+run_name: "NGP_20220327"
+
+# (Optional) a file of point list 
+# if you just want to run default pointing:
+# - pointing_dir: null
+# - pointing_file: null
+pointing_dir: null
+pointing_file: null
+
+# Whether to use MPI
+run_option:
+  use_mpi: YES
+  # NOTE: "n_threads" paramters is currently not used in the backend
+  # simulation codes. It should be implemented later in the web frontend
+  # in order to config the number of threads to request from NAOC cluster
+  n_threads: 80
+
+  # Output catalog only?
+  # If yes, no imaging simulation will run
+  out_cat_only: NO
+
+  # Only simulate stars?
+  star_only: NO
+
+  # Only simulate galaxies?
+  galaxy_only: NO
+
+###############################################
+# Observation setting
+###############################################
+obs_setting:
+
+  # Options for survey types:
+  # "Photometric": simulate photometric chips only
+  # "Spectroscopic": simulate slitless spectroscopic chips only
+  # "All": simulate full focal plane
+  survey_type: "All"
+  
+  # Exposure time [seconds]
+  exp_time: 150.
+
+  # Observation starting date & time
+  # (Subject to change)
+  date_obs: "210525" # [yymmdd]
+  time_obs: "120000" # [hhmmss]
+
+  # Default Pointing [degrees]
+  # Note: NOT valid when a pointing list file is specified
+  ra_center: 192.8595
+  dec_center:  27.1283
+
+  # Image rotation [degree]
+  image_rot: -113.4333
+
+  # Number of calibration pointings
+  np_cal: 0
+
+  # Run specific pointing(s):
+  # - give a list of indexes of pointings: [ip_1, ip_2...]
+  # - run all pointings: null
+  # Note: only valid when a pointing list is specified
+  run_pointings: null
+  
+  # Run specific chip(s):
+  # - give a list of indexes of chips: [ip_1, ip_2...]
+  # - run all chips: null
+  # Note: for all pointings
+  run_chips: null
+
+  # Whether to enable astrometric modeling
+  # astrometric_lib: "libshao.so"
+  enable_astrometric_model: True
+
+###############################################
+# Input path setting
+###############################################
+
+# Default path settings for WIDE survey simulation
+input_path:
+  cat_dir: "OnOrbitCalibration/CTargets20211231"
+  star_cat: "CT-NGP_V2.2.hdf5"
+  galaxy_cat: "galaxyCats_r_10.0_healpix_shift_192.859500_27.128300.hdf5"
+
+SED_templates_path:
+  star_SED: "Catalog_20210126/SpecLib.hdf5"
+  galaxy_SED: "Templates/Galaxy/"
+
+###############################################
+# PSF setting
+###############################################
+psf_setting:
+  
+  # Which PSF model to use:
+  # "Gauss": simple gaussian profile
+  # "Interp": Interpolated PSF from sampled ray-tracing data
+  psf_model: "Interp"
+
+  # PSF size [arcseconds]
+  # radius of 80% energy encircled
+  # NOTE: only valid for "Gauss" PSF
+  psf_rcont: 0.15
+
+  # path to PSF data
+  # NOTE: only valid for "Interp" PSF
+  psf_dir: "/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube"
+
+  # path to field-distortion model
+  # Note: only valid when ins_effects: field_dist is "ON"
+  fd_path: "FieldDistModelGlobal_v1.0.pickle"
+  
+  # sigma_spin: 0.0 # psf spin?
+
+###############################################
+# Shear setting
+###############################################
+
+shear_setting:
+  # Options to generate mock shear field:
+  # "constant": all galaxies are assigned a constant reduced shear
+  # "catalog": from catalog (not available yet)
+  # "extra": from seprate file
+  shear_type: "constant"
+
+  # For constant shear filed
+  reduced_g1: 0.026
+  reduced_g2: 0.015
+
+  # Representation of the shear vector?
+  reShear: "E"
+
+  # rotate galaxy ellipticity
+  rotateEll: 0. # [degree]
+
+  # Extra shear catalog
+  # (currently not used)
+  # shear_cat: "mockShear.cat"
+
+###############################################
+# Instrumental effects setting
+###############################################
+ins_effects:
+  # switches
+  field_dist:     ON  # Whether to add field distortions
+  add_back:       ON  # Whether to add sky background
+  add_dark:       ON  # Whether to add dark noise
+  add_readout:    ON # Whether to add read-out (Gaussian) noise
+  add_bias:       ON # Whether to add bias-level to images
+  shutter_effect: ON # Whether to add shutter effect
+  flat_fielding:  ON # Whether to add flat-fielding effect
+  prnu_effect:    ON # Whether to add PRNU effect
+  non_linear:     ON # Whether to add non-linearity
+  cosmic_ray:     ON # Whether to add cosmic-ray
+  cray_differ:    ON # Whether to generate different cosmic ray maps CAL and MS output
+  cte_trail:      OFF # Whether to simulate CTE trails
+  saturbloom:     ON # Whether to simulate Saturation & Blooming
+  add_badcolumns: ON # Whether to add bad columns
+  add_hotpixels:  ON # Whether to add hot pixels
+  add_deadpixels: ON # Whether to add dead(dark) pixels
+  bright_fatter:  ON # Whether to simulate Brighter-Fatter (also diffusion) effect
+
+  # values
+  dark_exptime:   300   # Exposure time for dark current frames [seconds]
+  flat_exptime:   150   # Exposure time for flat-fielding frames [seconds]
+  readout_time:   40    # The read-out time for each channel [seconds]
+  df_strength:    2.3   # Sillicon sensor diffusion strength
+  bias_level:     500   # bias level [e-/pixel]
+  gain:           1.1   # Gain
+  full_well:      90000 # Full well depth [e-]
+  NBias:          1     # Number of bias frames to be exported for each exposure
+  NDark:          1     # Number of dark frames to be exported for each exposure
+  NFlat:          1     # Number of flat frames to be exported for each exposure
+
+###############################################
+# Output options
+###############################################
+output_setting:
+  readout16:      OFF # Whether to export as 16 channels (subimages) with pre- and over-scan
+  shutter_output: OFF # Whether to export shutter effect 16-bit image
+  bias_output:    ON  # Whether to export bias frames
+  dark_output:    ON  # Whether to export the combined dark current files
+  flat_output:    ON  # Whether to export the combined flat-fielding files
+  prnu_output:    OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files
+
+###############################################
+# Random seeds
+###############################################
+random_seeds:
+  seed_Av:              121212    # Seed for generating random intrinsic extinction
+  seed_poisson:         20210601  # Seed for Poisson noise
+  seed_CR:              20210317  # Seed for generating random cosmic ray maps
+  seed_flat:            20210101  # Seed for generating random flat fields
+  seed_prnu:            20210102  # Seed for photo-response non-uniformity
+  seed_gainNonUniform:  20210202  # Seed for gain nonuniformity
+  seed_biasNonUniform:  20210203  # Seed for bias nonuniformity
+  seed_rnNonUniform:    20210204  # Seed for readout-noise nonuniformity
+  seed_badcolumns:      20240309  # Seed for bad columns
+  seed_defective:       20210304  # Seed for defective (bad) pixels
+  seed_readout:         20210601  # Seed for read-out gaussian noise
+...
\ No newline at end of file

run_NGP.pbs

+#!/bin/bash
+
+#PBS -N SIMS
+##PBS -l walltime=70:00:00
+
+##mpdallexit
+##mpdboot -n 10 -f  ./mpd.hosts
+
+##PBS -j oe
+#PBS -l nodes=comput110:ppn=80
+#####PBS -q longq
+#PBS -q batch
+#PBS -u fangyuedong
+
+NP=30
+date
+echo $NP
+
+mpirun -np $NP python /public/home/fangyuedong/temp/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/temp/CSST/config

tests/PSFInterpTest/PSFInterpModule_coverage.py

+import unittest
+
+import sys,os,math
+from itertools import islice
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+mpl.use('Agg')
+
+import yaml
+from ObservationSim.Config import Config
+from ObservationSim.Config.Config import config_dir
+from ObservationSim.Instrument import Chip
+
+from ObservationSim.PSF.PSFInterp import PSFInterp
+
+
+def defineCCD(iccd, config_file):
+    with open(config_file, "r") as stream:
+        try:
+            config = yaml.safe_load(stream)
+            #for key, value in config.items():
+            #    print (key + " : " + str(value))
+        except yaml.YAMLError as exc:
+            print(exc)
+    path_dict = config_dir(config=config, work_dir=config['work_dir'], data_dir=config['data_dir'])
+    chip = Chip(chipID=iccd, config=config)
+    #chip = Chip(chipID=iccd, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict['sls_dir'], config=config)
+    return chip
+
+def psfSecondMoments(psfMat, cenX, cenY, pixSize=1):
+    apr = 0.5 #arcsec, 0.5角秒内测量
+    fl  = 28. #meters
+    pxs = 2.5 #microns
+    apr = np.deg2rad(apr/3600.)*fl*1e6
+    apr = apr/pxs
+    apr = np.int(np.ceil(apr))
+
+    I = psfMat
+    ncol = I.shape[1]
+    nrow = I.shape[0]
+    w   = 0.0
+    w11 = 0.0
+    w12 = 0.0
+    w22 = 0.0
+    for icol in range(ncol):
+        for jrow in range(nrow):
+            x = icol*pixSize - cenX
+            y = jrow*pixSize - cenY
+            rr = np.sqrt(x*x + y*y)
+            wgt= 0.0
+            if rr <= apr:
+                wgt = 1.0
+            w   += I[jrow, icol]*wgt
+            w11 += x*x*I[jrow, icol]*wgt
+            w12 += x*y*I[jrow, icol]*wgt
+            w22 += y*y*I[jrow, icol]*wgt
+    w11 /= w
+    w12 /= w
+    w22 /= w
+    sz = w11 + w22
+    e1 = (w11 - w22)/sz
+    e2 = 2.0*w12/sz
+
+    return sz, e1, e2
+
+
+def test_psfEll(iccd, iwave, psfMat):
+    psfMat_iwave = psfMat.psfMat[iwave-1, :,:,:]
+    npsf = np.shape(psfMat_iwave)[0]
+
+    imx = np.zeros(npsf)
+    imy = np.zeros(npsf)
+    psf_e1 = np.zeros(npsf)
+    psf_e2 = np.zeros(npsf)
+    psf_sz = np.zeros(npsf)
+
+    for ipsf in range(1, npsf+1):
+        print('ipsf-{:}'.format(ipsf), end='\r')
+        imx[ipsf-1] = psfMat.cen_col[iwave-1, ipsf-1]
+        imy[ipsf-1] = psfMat.cen_row[iwave-1, ipsf-1]
+        psfMat_iwave_ipsf = psfMat_iwave[ipsf-1, :, :]
+        cenX = 256
+        cenY = 256
+        sz, e1, e2  = psfSecondMoments(psfMat_iwave_ipsf, cenX, cenY, pixSize=1)
+        psf_e1[ipsf-1] = e1
+        psf_e2[ipsf-1] = e2
+        psf_sz[ipsf-1] = sz
+        #print('ell======', ipsf, np.sqrt(e1**2 + e2**2))
+    #######
+    arr = [imx, imy, psf_e1, psf_e2, psf_sz]
+    np.save('data/psfEll{:}_{:}_{:}'.format(int(np.sqrt(npsf)),iccd, iwave), arr)
+
+
+def test_psfEllPlot(OVERPLOT=False):
+    #if ThisTask == 0:
+    if True:
+        prefix = 'psfEll30'
+        iccd = 1
+        iwave= 1
+        data = np.load('data/'+prefix+'_1_1.npy')
+        imx= data[0]
+        imy= data[1]
+        psf_e1 = data[2]
+        psf_e2 = data[3]
+        print(np.shape(imx))
+        npsf = np.shape(imx)[0]
+
+        #######
+        plt.cla()
+        plt.close("all")
+
+        fig = plt.figure(figsize=(12, 12))
+        plt.subplots_adjust(wspace=0.1, hspace=0.1)
+        ax = plt.subplot(1, 1, 1)
+        for ipsf in range(npsf):
+            plt.plot(imx[ipsf], imy[ipsf], 'r.')
+
+            ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2
+            ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2)
+            ell *= 15
+            lcos = ell*np.cos(ang)
+            lsin = ell*np.sin(ang)
+            plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=2)
+            ###########
+        ang = 0.
+        ell = 0.05
+        ell*= 15
+        lcos = ell*np.cos(ang)
+        lsin = ell*np.sin(ang)
+        plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2)
+        plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10)
+        plt.xlabel('CCD X (mm)')
+        plt.ylabel('CCD Y (mm)')
+
+        if OVERPLOT == True:
+            prefix = 'psfEll20'
+            data = np.load('data/'+prefix+'_1_1.npy')
+            imx= data[0]
+            imy= data[1]
+            psf_e1 = data[2]
+            psf_e2 = data[3]
+
+            npsf = np.shape(imx)[0]
+            for ipsf in range(npsf):
+                plt.plot(imx[ipsf], imy[ipsf], 'b.')
+
+                ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2
+                ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2)
+                ell *= 15
+                lcos = ell*np.cos(ang)
+                lsin = ell*np.sin(ang)
+                plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2)
+
+        plt.gca().set_aspect(1)
+        if OVERPLOT == True:
+            prefix = 'psfEllOP'
+        plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd))
+
+
+def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB):
+    bandpass = iwave-1
+    class pos_img():
+        def __init__(self,x, y):
+            self.x = x*1e3/10.  #in unit of pixels
+            self.y = y*1e3/10.
+
+    psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:]
+    npsf = np.shape(psfMat_iwave)[0]
+    psf_e1 = np.zeros(npsf)
+    psf_e2 = np.zeros(npsf)
+    psf_sz = np.zeros(npsf)
+
+    for ipsf in range(1, npsf+1):
+        print('ipsf:', ipsf, end='\r', flush=True)
+        tpos_img = pos_img(psfMatA.cen_col[iwave-1, ipsf-1], psfMatA.cen_row[iwave-1, ipsf-1])
+        psfIDW = psfMatB.get_PSF(chip, tpos_img, bandpass, galsimGSObject=False, findNeighMode='treeFind')
+        np.save('figs/psfIDW/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf), psfIDW)
+
+        cenX = 256
+        cenY = 256
+        sz, e1, e2  = psfSecondMoments(psfIDW, cenX, cenY, pixSize=1)
+        psf_e1[ipsf-1] = e1
+        psf_e2[ipsf-1] = e2
+        psf_sz[ipsf-1] = sz
+    arr = [psf_e1, psf_e2, psf_sz]
+    np.save('data/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr)
+
+
+def test_psfResidualPlot(iccd, iwave, ipsf, psfMatA):
+    psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:]
+
+    psfMatORG = psfMat_iwave[ipsf-1, :, :]
+    psfMatIDW = np.load('figs/psfIDW/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf))
+
+    npix = psfMatORG.shape[0]
+    pixCutEdge= int(npix/2-15)
+
+    img0 = psfMatORG[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge]
+    img1 = psfMatIDW[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge]
+    imgX = (img1 - img0)/img0
+
+    img0 = np.log10(img0)
+    img1 = np.log10(img1)
+    imgX = np.log10(np.abs(imgX))
+
+    fig = plt.figure(figsize=(18,4))
+    ax = plt.subplot(1,3,1)
+    plt.imshow(img0, origin='lower', vmin=-7, vmax=-1.3)
+    plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--')
+    plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--')
+    plt.annotate('ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15)
+    cticks=[-7, -6, -5, -4, -3, -2, -1]
+    cbar = plt.colorbar(ticks=cticks)
+    cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$'])
+    print(img0.min(), img0.max())
+
+
+    ax = plt.subplot(1,3,2)
+    plt.imshow(img1, origin='lower', vmin=-7, vmax=-1.3)
+    plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--')
+    plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--')
+    plt.annotate('IDW', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15)
+    cticks=[-7, -6, -5, -4, -3, -2, -1]
+    cbar = plt.colorbar(ticks=cticks)
+    cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$'])
+    print(img1.min(), img1.max())
+
+
+    ax = plt.subplot(1,3,3)
+    plt.imshow(imgX, origin='lower', vmin =-3, vmax =np.log10(3e-1))
+    plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--')
+    plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--')
+    #plt.annotate('(IDW-ORG)/ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15)
+    cticks=[-5, -4, -3, -2, -1]
+    cbar = plt.colorbar(ticks=cticks)
+    cbar.ax.set_yticklabels(['$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$'])
+
+    print(np.max((psfMatORG-psfMatIDW)))
+    plt.savefig('figs/psfResidual_iccd{:}.pdf'.format(iccd))
+
+
+def test_psfEllIDWPlot(OVERPLOT=False):
+    #if ThisTask == 0:
+    if True:
+        prefix = 'psfEll20'
+        iccd = 1
+        iwave= 1
+        data = np.load('data/'+prefix+'_1_1.npy')
+        imx= data[0]
+        imy= data[1]
+        psf_e1 = data[2]
+        psf_e2 = data[3]
+        print(np.shape(imx))
+        npsf = np.shape(imx)[0]
+
+        #######
+        plt.cla()
+        plt.close("all")
+
+        fig = plt.figure(figsize=(12, 12))
+        plt.subplots_adjust(wspace=0.1, hspace=0.1)
+        ax = plt.subplot(1, 1, 1)
+        for ipsf in range(npsf):
+            plt.plot(imx[ipsf], imy[ipsf], 'b.')
+
+            ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2
+            ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2)
+            ell *= 15
+            lcos = ell*np.cos(ang)
+            lsin = ell*np.sin(ang)
+            plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2)
+            ###########
+        ang = 0.
+        ell = 0.05
+        ell*= 15
+        lcos = ell*np.cos(ang)
+        lsin = ell*np.sin(ang)
+        #plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2)
+        #plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10)
+        plt.xlabel('CCD X (mm)')
+        plt.ylabel('CCD Y (mm)')
+
+        if OVERPLOT == True:
+            prefix = 'psfEll20IDW'
+            data = np.load('data/'+prefix+'_1_1.npy')
+            #imx= data[0]
+            #imy= data[1]
+            psf_e1 = data[0]
+            psf_e2 = data[1]
+
+            npsf = np.shape(imx)[0]
+            for ipsf in range(npsf):
+                #plt.plot(imx[ipsf], imy[ipsf], 'r.')
+
+                ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2
+                ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2)
+                ell *= 15
+                lcos = ell*np.cos(ang)
+                lsin = ell*np.sin(ang)
+                plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=1)
+
+        plt.gca().set_aspect(1)
+        if OVERPLOT == True:
+            prefix = 'psfEllOPIDW'
+        plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd))
+
+
+def test_psfdEllabsPlot(iccd):
+    #if ThisTask == 0:
+    if True:
+        prefix = 'psfEll20'
+        #iccd = 1
+        #iwave= 1
+        data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd))
+        imx= data[0]
+        imy= data[1]
+        psf_e1 = data[2]
+        psf_e2 = data[3]
+        psf_sz = data[4]
+        print(np.shape(imx))
+        npsf = np.shape(imx)[0]
+
+        ellX = np.sqrt(psf_e1**2 + psf_e2**2)
+        angX = np.arctan2(psf_e2, psf_e1)/2
+        angX = np.rad2deg(angX)
+        szX  = psf_sz
+
+        ##############################
+        prefix = 'psfEll20IDW'
+        data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd))
+        #imx= data[0]
+        #imy= data[1]
+        psf_e1 = data[0]
+        psf_e2 = data[1]
+        psf_sz = data[2]
+
+        ellY = np.sqrt(psf_e1**2 + psf_e2**2)
+        angY = np.arctan2(psf_e2, psf_e1)/2
+        angY = np.rad2deg(angY)
+        szY  = psf_sz
+
+        ##############################
+        fig=plt.figure(figsize=(6, 5))
+        grid = plt.GridSpec(3,1,left=0.15, right=0.95, wspace=None, hspace=0.02)
+        #plt.subplots_adjust(left=None,bottom=None,right=None,top=None,wspace=None,hspace=0.02)
+        ax = plt.subplot(grid[0:2,0])
+        plt.plot([0.01,0.1],[0.01,0.1], 'k--', lw=1. )
+        plt.scatter(ellX, ellY, color='b', alpha=1., s=3., edgecolors='None')
+        #plt.xlim([0.015, 0.085])
+        #plt.ylim([0.015, 0.085])
+        plt.ylabel('$\epsilon_{\\rm IDW}$')
+        plt.gca().axes.get_xaxis().set_visible(False)
+
+        ax = plt.subplot(grid[2,0])
+        plt.plot([0.015,0.085],[0.,0.], 'k--', lw=1. )
+        plt.scatter(ellX, (ellY-ellX), color='b', s=3., edgecolors='None')
+        #plt.xlim([0.015, 0.085])
+        #plt.ylim([-0.0018, 0.0018])
+        plt.xlabel('$\epsilon_{\\rm ORG}$')
+        plt.ylabel('$\Delta$')
+        plt.savefig('figs/psfEllOPIDWPDF_{:}.pdf'.format(iccd))
+
+        fig=plt.figure(figsize=(6, 6))
+        plt.hist((szY-szX)/szX, bins=20, color='r', alpha=0.5)
+        plt.xlabel('$(R_{\\rm IDW}-R_{\\rm ORG})/R_{\\rm ORG}$')
+        plt.ylabel('PDF')
+        plt.savefig('figs/psfEllOPIDWPDF_dsz_{:}.pdf'.format(iccd))
+
+
+class PSFInterpModule_coverage(unittest.TestCase):
+    def test_psfEll_(self):
+        iccd = 1
+        iwave= 1
+
+        config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml"
+        chip = defineCCD(iccd, config_file)
+        print(chip.chipID)
+        print(chip.cen_pix_x, chip.cen_pix_y)
+
+        ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCubeTest'
+        psfMatA = PSFInterp(chip, npsf=400, PSF_data_file=ipath, PSF_data_prefix="S20x20_")
+        psfMatB = PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="S30x30_")
+
+        test_psfEll(iccd, iwave, psfMatA)
+        test_psfEll(iccd, iwave, psfMatB)
+        test_psfEllPlot(OVERPLOT=True)
+
+        test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB)
+        ipsf = 1
+        test_psfResidualPlot(iccd, iwave, ipsf, psfMatA)
+        test_psfEllIDWPlot(OVERPLOT=True)
+        test_psfdEllabsPlot(iccd)
+
+
+if __name__ == '__main__':
+    unittest.main()
+    print('#####haha#####')
+
+
+

tests/PSFInterpTest/loadPSFSet.py

+import unittest
+
+import sys,os,math
+from itertools import islice
+
+import numpy as np
+
+import yaml
+from ObservationSim.Config import Config
+from ObservationSim.Config.Config import config_dir
+from ObservationSim.Instrument import Chip
+
+from ObservationSim.PSF.PSFInterp import PSFInterp
+
+
+def defineCCD(iccd, config_file):
+    with open(config_file, "r") as stream:
+        try:
+            config = yaml.safe_load(stream)
+            #for key, value in config.items():
+            #    print (key + " : " + str(value))
+        except yaml.YAMLError as exc:
+            print(exc)
+    path_dict = config_dir(config=config, work_dir=config['work_dir'], data_dir=config['data_dir'])
+    chip = Chip(chipID=iccd, config=config)
+    #chip = Chip(chipID=iccd, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict['sls_dir'], config=config)
+    return chip
+
+def loadPSFSet(iccd):
+    config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml"
+    chip = defineCCD(iccd, config_file)
+    print(chip.chipID)
+    print(chip.cen_pix_x, chip.cen_pix_y)
+
+    ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube'
+    psfMat= PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="")
+    psfSet= psfMat._loadPSF(iccd, ipath, PSF_data_prefix="")
+
+    twave = 0 #[0...3]
+    tpsf  = 0 #[0...899]
+    field_x= psfSet[twave][tpsf]['field_x']
+    field_y= psfSet[twave][tpsf]['field_y']
+    image_x= psfSet[twave][tpsf]['image_x']
+    image_y= psfSet[twave][tpsf]['image_y']
+    centroid_x= psfSet[twave][tpsf]['centroid_x']
+    centroid_y= psfSet[twave][tpsf]['centroid_y']
+    print("pos_info:", field_x, field_y, image_x, image_y, centroid_x, centroid_y)
+    return psfSet
+
+
+
+class PSFInterpModule_coverage(unittest.TestCase):
+    def test_psfEll_(self):
+        iccd = 1  #[1...30]
+        psfSet = loadPSFSet(iccd)
+
+
+if __name__ == '__main__':
+    unittest.main()
+    print('#####haha#####')
+
+
+