Merge branch 'wcs_test' into 'develop'

starting point of the new version

See merge request csst_sim/csst-simulation!15

Catalog/C6_50sqdeg.py

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

Catalog/C6_Catalog.py

@@ -83,7 +83,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -259,7 +259,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/C6_fits_Catalog.py

@@ -94,7 +94,7 @@ class Catalog(CatalogBase):
         ###mock_stamp_END
         
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -272,7 +272,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/FGS_Catalog.py

@@ -84,7 +84,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -260,7 +260,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/NGPCatalog.py

@@ -52,7 +52,7 @@ class Catalog(CatalogBase):
             self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
             self._load_SED_lib_gals()
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -191,7 +191,7 @@ class Catalog(CatalogBase):
             param['id'] = gals['galaxyID'][igals]
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             if param['star'] == 2:
                 obj = Quasar(param, logger=self.logger)
             

Catalog/fd_test_C6.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+from astropy.coordinates import spherical_to_cartesian
+from astropy.table import Table
+from scipy import interpolate
+from datetime import datetime
+
+from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=ngals,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        # for igals in range(ngals):
+        for igals in range(0, ngals, 5):
+            # # (TEST)
+            # if igals > 1000:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            param['mag_use_normal'] = 20.
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            # param['g1'] = gals['shear'][igals][0]
+            # param['g2'] = gals['shear'][igals][1]
+            param['g1'] = 0.
+            param['g2'] = 0.
+            param['kappa'] = gals['kappa'][igals]
+            param['e1'] = gals['ellipticity_true'][igals][0]
+            param['e2'] = gals['ellipticity_true'][igals][1]
+            
+            # For shape calculation
+            
+            param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            # param['e1_disk'] = param['e1']
+            # param['e2_disk'] = param['e2']
+            # param['e1_bulge'] = param['e1']
+            # param['e2_bulge'] = param['e2']
+            param['e1_disk'] = 0.
+            param['e2_disk'] = 0.
+            param['e1_bulge'] = 0.
+            param['e2_bulge'] = 0.
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            # param['size'] = gals['size'][igals]
+            param['size'] = 1.
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            # param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            param['bfrac'] = 0.
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["Dec"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=nstars,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_arr[istars]
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            param['mag_use_normal'] = stars['app_sdss_g'][istars]
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=nAGNs,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        # if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+        #     star_cat = h5.File(self.star_path, 'r')['catalog']
+        #     for pix in self.pix_list:
+        #         try:
+        #             stars = star_cat[str(pix)]
+        #             self._load_stars(stars, pix_id=pix)
+        #             del stars
+        #         except Exception as e:
+        #             self.logger.error(str(e))
+        #             print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        # if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+        #     try:
+        #         self._load_AGNs()
+        #     except Exception as e:
+        #         traceback.print_exc()
+        #         self.logger.error(str(e))
+        #         print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/testCat_galaxy.py

@@ -81,7 +81,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -253,7 +253,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/wcs_test_C6.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+from astropy.coordinates import spherical_to_cartesian
+from astropy.table import Table
+from scipy import interpolate
+from datetime import datetime
+
+from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
+from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=ngals,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        for igals in range(ngals):
+            # # (TEST)
+            # if igals > 100:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            # param['mag_use_normal'] = 20.
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            # param['g1'] = gals['shear'][igals][0]
+            # param['g2'] = gals['shear'][igals][1]
+            param['g1'] = 0.
+            param['g2'] = 0.
+            param['kappa'] = gals['kappa'][igals]
+            # param['e1'] = gals['ellipticity_true'][igals][0]
+            # param['e2'] = gals['ellipticity_true'][igals][1]
+            param['e1'] = 0.
+            param['e2'] = 0.
+            
+            # For shape calculation
+            
+            param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            param['e1_disk'] = param['e1']
+            param['e2_disk'] = param['e2']
+            param['e1_bulge'] = param['e1']
+            param['e2_bulge'] = param['e2']
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            param['size'] = gals['size'][igals]
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["Dec"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=nstars,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_arr[istars]
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            # param['mag_use_normal'] = stars['app_sdss_g'][istars]
+            param['mag_use_normal'] = 20.
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            dt = datetime.utcfromtimestamp(self.pointing.timestamp)
+            date_str = dt.date().isoformat()
+            time_str = dt.time().isoformat()
+            ra_arr, dec_arr = on_orbit_obs_position(
+                input_ra_list=ra_list,
+                input_dec_list=dec_list,
+                input_pmra_list=pmra_list,
+                input_pmdec_list=pmdec_list,
+                input_rv_list=rv_list,
+                input_parallax_list=parallax_list,
+                input_nstars=nAGNs,
+                input_x=self.pointing.sat_x,
+                input_y=self.pointing.sat_y,
+                input_z=self.pointing.sat_z,
+                input_vx=self.pointing.sat_vx,
+                input_vy=self.pointing.sat_vy,
+                input_vz=self.pointing.sat_vz,
+                input_epoch="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['catalog']
+            for pix in self.pix_list:
+                try:
+                    stars = star_cat[str(pix)]
+                    self._load_stars(stars, pix_id=pix)
+                    del stars
+                except Exception as e:
+                    self.logger.error(str(e))
+                    print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+            try:
+                self._load_AGNs()
+            except Exception as e:
+                traceback.print_exc()
+                self.logger.error(str(e))
+                print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

ObservationSim/Config/Header/ImageHeader.py

@@ -291,7 +291,8 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
         r_dat['CRPIX1'] = -x1
         r_dat['CRPIX2'] = -y1
 
-        cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
+        # cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
+        cd = np.array([[ pixel_scale,  0], [0, -pixel_scale]])/3600.
         cd_rot = rotate_CD_matrix(cd, pa_aper)
         r_dat['CD1_1'] = cd_rot[0,0]
         r_dat['CD1_2'] = cd_rot[0,1]

ObservationSim/Instrument/Chip/Chip.py

@@ -77,7 +77,8 @@ class Chip(FocalPlane):
                 self.fdModel = None
             else:
                 try:
-                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                    # with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion:
                         with open(field_distortion, "rb") as f:
                             self.fdModel = pickle.load(f)
                 except AttributeError:
@@ -117,7 +118,7 @@ class Chip(FocalPlane):
         self._getCRdata()
 
         # Define the sensor model
-        if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
+        if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
             self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
         else:
             self.sensor = galsim.Sensor()
@@ -146,18 +147,18 @@ class Chip(FocalPlane):
     def _getSurveyType(self):
         if self.filter_type in ["GI", "GV", "GU"]:
             return "spectroscopic"
-        elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
+        elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
             return "photometric"
         # elif self.filter_type in ["FGS"]:
         #     return "FGS"
 
     def _getChipEffCurve(self, filter_type):
         # CCD efficiency curves
-        if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
+        if filter_type in ['NUV', 'u', 'GU']: filename = 'UV0.txt'
         if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
         if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
         # Mirror efficiency:
-        # if filter_type == 'nuv': mirror_eff = 0.54
+        # if filter_type == 'NUV': mirror_eff = 0.54
         # if filter_type == 'u': mirror_eff = 0.68
         # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
         # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
@@ -184,7 +185,7 @@ class Chip(FocalPlane):
     def getChipFilter(self, chipID=None, filter_layout=None):
         """Return the filter index and type for a given chip #(chipID)
         """
-        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
+        filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
         if filter_layout is not None:
             return filter_layout[chipID][0], filter_layout[chipID][1]
         if chipID == None:
@@ -197,7 +198,7 @@ class Chip(FocalPlane):
         if chipID in [7, 24]:          filter_type = "i"
         if chipID in [14, 17]:         filter_type = "u"
         if chipID in [9, 22]:          filter_type = "r"
-        if chipID in [12, 13, 18, 19]: filter_type = "nuv"
+        if chipID in [12, 13, 18, 19]: filter_type = "NUV"
         if chipID in [8, 23]:          filter_type = "g"
         if chipID in [1, 10, 21, 30]:  filter_type = "GI"
         if chipID in [2, 5, 26, 29]:   filter_type = "GV"
@@ -514,21 +515,6 @@ class Chip(FocalPlane):
         if config["ins_effects"]["add_badcolumns"] == True:
             img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
 
-        # Add Bias level
-        if config["ins_effects"]["add_bias"] == True:
-            if self.logger is not None:
-                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
-            else:
-                print("  Adding Bias level and 16-channel non-uniformity")
-            if config["ins_effects"]["bias_16channel"] == True:
-                img = effects.AddBiasNonUniform16(img, 
-                    bias_level=float(self.bias_level), 
-                    nsecy = 2, nsecx=8, 
-                    seed=SeedBiasNonuni+self.chipID,
-                    logger=self.logger)
-            elif config["ins_effects"]["bias_16channel"] == False:
-                img += self.bias_level
-
         # Apply Nonlinearity on the chip image
         if config["ins_effects"]["non_linear"] == True:
             if self.logger is not None:
@@ -553,6 +539,21 @@ class Chip(FocalPlane):
                 print("  Apply CTE Effect")
             img = effects.CTE_Effect(GSImage=img, threshold=27)
         
+        # Add Bias level
+        if config["ins_effects"]["add_bias"] == True:
+            if self.logger is not None:
+                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
+            else:
+                print("  Adding Bias level and 16-channel non-uniformity")
+            if config["ins_effects"]["bias_16channel"] == True:
+                img = effects.AddBiasNonUniform16(img, 
+                    bias_level=float(self.bias_level), 
+                    nsecy = 2, nsecx=8, 
+                    seed=SeedBiasNonuni+self.chipID,
+                    logger=self.logger)
+            elif config["ins_effects"]["bias_16channel"] == False:
+                img += self.bias_level
+
         # Add Read-out Noise
         if config["ins_effects"]["add_readout"] == True:
             seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID

ObservationSim/Instrument/FilterParam.py

@@ -26,7 +26,7 @@ class FilterParam(object):
         # 8) dim end magnitude
         if filter_param == None:
             filtP = {
-                    "nuv":   [2867.7,   705.4,  2470.0,   3270.0,  0.1404, 0.004, 15.7, 25.4],
+                    "NUV":   [2867.7,   705.4,  2470.0,   3270.0,  0.1404, 0.004, 15.7, 25.4],
                     "u":     [3601.1,   852.1,  3120.0,   4090.0,  0.2176, 0.021, 16.1, 25.4],
                     "g":     [4754.5,  1569.8,  3900.0,   5620.0,  0.4640, 0.164, 17.2, 26.3],
                     "r":     [6199.8,  1481.2,  5370.0,   7030.0,  0.5040, 0.207, 17.0, 26.0],

ObservationSim/Instrument/FocalPlane.py

@@ -86,15 +86,16 @@ class FocalPlane(object):
         if (xcen == None) or (ycen == None):
             xcen = self.cen_pix_x
             ycen = self.cen_pix_y
+        # dudx =  -np.cos(img_rot.rad) * pix_scale
+        # dudy =  -np.sin(img_rot.rad) * pix_scale
+        # dvdx =  -np.sin(img_rot.rad) * pix_scale
+        # dvdy =  +np.cos(img_rot.rad) * pix_scale
+
         dudx =  -np.cos(img_rot.rad) * pix_scale
-        dudy =  -np.sin(img_rot.rad) * pix_scale
+        dudy =  +np.sin(img_rot.rad) * pix_scale
         dvdx =  -np.sin(img_rot.rad) * pix_scale
-        dvdy =  +np.cos(img_rot.rad) * pix_scale
+        dvdy =  -np.cos(img_rot.rad) * pix_scale
         
-        # dudx =  +np.sin(img_rot.rad) * pix_scale
-        # dudy =  +np.cos(img_rot.rad) * pix_scale
-        # dvdx =  -np.cos(img_rot.rad) * pix_scale
-        # dvdy =  +np.sin(img_rot.rad) * pix_scale
         moscen = galsim.PositionD(x=xcen, y=ycen)
         sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)
         affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)

ObservationSim/MockObject/CatalogBase.py

 import numpy as np
 import galsim
 import copy
+import cmath
 
 from astropy.table import Table
 from abc import abstractmethod, ABCMeta
@@ -73,6 +74,16 @@ class CatalogBase(metaclass=ABCMeta):
         }
         return param
 
+    @staticmethod
+    def rotate_ellipticity(e1, e2, rotation=0., unit='radians'):
+        if unit == 'degree':
+            rotation = np.radians(rotation)
+        e_total = np.sqrt(e1**2 + e2**2)
+        phi = cmath.phase(complex(e1, e2))
+        e1 = e_total * np.cos(phi + 2*rotation)
+        e2 = e_total * np.sin(phi + 2*rotation)
+        return e1, e2, e_total
+    
     @staticmethod
     def convert_sed(mag, sed, target_filt, norm_filt=None):
         bandpass = target_filt.bandpass_full

ObservationSim/MockObject/Galaxy.py

@@ -12,7 +12,7 @@ from ObservationSim.MockObject.MockObject import MockObject
 # import tracemalloc
 
 class Galaxy(MockObject):
-    def __init__(self, param, rotation=None, logger=None):
+    def __init__(self, param, logger=None):
         super().__init__(param, logger=logger)
         # self.thetaR = self.param["theta"]
         # self.bfrac = self.param["bfrac"]
@@ -27,8 +27,6 @@ class Galaxy(MockObject):
         # 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)
         if not hasattr(self, "disk_sersic_idx"):
             self.disk_sersic_idx = 1.
         if not hasattr(self, "bulge_sersic_idx"):
@@ -51,6 +49,7 @@ class Galaxy(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            if self.logger:
                 self.logger.error(e)
             return -1
         for i in range(len(bandpass_list)):
@@ -59,6 +58,7 @@ class Galaxy(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                if self.logger:
                     self.logger.error(e)
                 return -1
             
@@ -78,12 +78,13 @@ class Galaxy(MockObject):
 
             gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
             gal = gal.withFlux(nphotons)
+            if fd_shear is not None:
+                g1 += fd_shear.g1
+                g2 += fd_shear.g2
             gal_shear = galsim.Shear(g1=g1, g2=g2)
             gal = gal.shear(gal_shear)
             gal = galsim.Convolve(psf, gal)
             
-            if fd_shear is not None:
-                gal = gal.shear(fd_shear)
             objs.append(gal)
         final = galsim.Sum(objs)
         return final
@@ -97,6 +98,7 @@ class Galaxy(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            if self.logger:
                 self.logger.error(e)
             return 2, None
 
@@ -131,6 +133,18 @@ class Galaxy(MockObject):
 
         real_wcs_local = self.real_wcs.local(self.real_pos)
 
+        disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
+        disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
+        disk = disk.shear(disk_shape)
+        bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
+        bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
+        bulge = bulge.shear(bulge_shape)
+        
+        if fd_shear:
+            g1 += fd_shear.g1
+            g2 += fd_shear.g2
+        gal_shear = galsim.Shear(g1=g1, g2=g2)
+
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
 
@@ -138,6 +152,7 @@ class Galaxy(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                if self.logger:
                     self.logger.error(e)
                 # return False
                 continue
@@ -153,41 +168,24 @@ class Galaxy(MockObject):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
-            disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
-            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
-            disk = disk.shear(disk_shape)
-            bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
-            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_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
+            gal_temp = gal_temp.shear(gal_shear)
+
+            gal_temp = gal_temp.withFlux(nphotons)
+            if not big_galaxy: # Not apply PSF for very big galaxy
+                gal_temp = galsim.Convolve(psf, gal_temp)
+            
+            if i == 0:
+                gal = gal_temp
+            else:
+                gal = gal + gal_temp
 
             # (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)
-            gal = gal.shear(gal_shear)
-
-            if not big_galaxy: # Not apply PSF for very big galaxy
-                gal = galsim.Convolve(psf, gal)
-                if fd_shear is not None:
-                    gal = gal.shear(fd_shear)
-
-            # Use (explicit) stamps to draw
-            stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
-            
-            xmax = max(xmax, stamp.xmax - stamp.xmin)
-            ymax = max(ymax, stamp.ymax - stamp.ymin)
-
-            photons = stamp.photons
-            photons.x += x_nominal
-            photons.y += y_nominal
-            photons_list.append(photons)
-            del gal
-
         # # [C6 TEST]
         # print('xmax = %d, ymax = %d '%(xmax, ymax))
         # # Output memory usage
@@ -196,7 +194,12 @@ class Galaxy(MockObject):
         # for stat in top_stats[:10]:
         #     print(stat)
 
-        stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1))
+        stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
+        photons = stamp.photons
+        photons.x += x_nominal
+        photons.y += y_nominal
+        photons_list.append(photons)
+
         stamp.wcs = real_wcs_local
         stamp.setCenter(x_nominal, y_nominal)
         bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
@@ -226,6 +229,7 @@ class Galaxy(MockObject):
         else:
             # Return code 0: object photons missed this detector
             print("obj %s missed"%(self.id))
+            if self.logger:
                 self.logger.info("obj %s missed"%(self.id))
             return 0, pos_shear
         
@@ -297,14 +301,17 @@ class Galaxy(MockObject):
             # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
 
             gal = gal.withFlux(tel.pupil_area * exptime)
+            if fd_shear:
+                g1 += fd_shear.g1
+                g2 += fd_shear.g2
             gal_shear = galsim.Shear(g1=g1, g2=g2)
             gal = gal.shear(gal_shear)
             gal = galsim.Convolve(psf, gal)
 
             if not big_galaxy: # Not apply PSF for very big galaxy
                 gal = galsim.Convolve(psf, gal)
-                if fd_shear is not None:
-                    gal = gal.shear(fd_shear)
+                # if fd_shear is not None:
+                #     gal = gal.shear(fd_shear)
 
             starImg = gal.drawImage(wcs=real_wcs_local, offset=offset)
 
@@ -397,14 +404,6 @@ class Galaxy(MockObject):
         final = galsim.Convolve(psf, gal)
         return final
 
-    def rotateEllipticity(self, rotation):
-        if rotation == 1:
-            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e2_disk, self.e1_disk, -self.e2_bulge, self.e1_bulge, -self.e2_total, self.e1_total
-        if rotation == 2:
-            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e1_disk, -self.e2_disk, -self.e1_bulge, -self.e2_bulge, -self.e1_total, -self.e2_total
-        if rotation == 3:
-            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = self.e2_disk, -self.e1_disk, self.e2_bulge, -self.e1_bulge, self.e2_total, -self.e1_total
-
     def drawObject(self, img, final, noise_level=0.0, flux=None, filt=None, tel=None, exptime=150.):
         """ Override the method in parent class 
         Need to constrain the size of image stamp for extended objects

ObservationSim/MockObject/MockObject.py

@@ -32,6 +32,7 @@ class MockObject(object):
 
         # Place holder for outputs
         self.additional_output_str = ""
+        
         self.fd_shear = None
 
         self.logger = logger
@@ -61,7 +62,7 @@ class MockObject(object):
         dec = self.param["dec"]
         return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees)
 
-    def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, img_header=None):
+    def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None):
         self.posImg = img.wcs.toImage(self.getPosWorld())
         self.localWCS = img.wcs.local(self.posImg)
         if (fdmodel is not None) and (chip is not None):
@@ -80,7 +81,9 @@ class MockObject(object):
         dy = y - self.y_nominal
         self.offset = galsim.PositionD(dx, dy)
         
-        if img_header is not None:
+        if chip_wcs is not None:
+            self.real_wcs = chip_wcs
+        elif img_header is not None:
             self.real_wcs = galsim.FitsWCS(header=img_header)
         else:
             self.real_wcs = None
@@ -132,6 +135,7 @@ class MockObject(object):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            if self.logger:
                 self.logger.error(e)
             return 2, None
 
@@ -164,6 +168,7 @@ class MockObject(object):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                if self.logger:
                     self.logger.error(e)
                 # return False
                 continue
@@ -214,6 +219,7 @@ class MockObject(object):
         else:
             # Return code 0: object photons missed this detector
             print("obj %s missed"%(self.id))
+            if self.logger:
                 self.logger.info("obj %s missed"%(self.id))
             return 0, pos_shear