diff --git a/Catalog/C6_50sqdeg.py b/Catalog/C6_50sqdeg.py new file mode 100644 index 0000000000000000000000000000000000000000..f13faa7f2df61a4205e49ef6d0def8301ef196a3 --- /dev/null +++ b/Catalog/C6_50sqdeg.py @@ -0,0 +1,470 @@ +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 diff --git a/Catalog/C6_Catalog.py b/Catalog/C6_Catalog.py index 5cef969ebff3223af77139a524c39641a0a9b6b3..7d114e573747b84a25e892e4938d02e020e355c8 100644 --- a/Catalog/C6_Catalog.py +++ b/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., diff --git a/Catalog/C6_fits_Catalog.py b/Catalog/C6_fits_Catalog.py index 47160dfe018ac19b1e008769dc0900c474e4b848..80af1cbe05224ff92a71424c4e6d859e6f957350 100644 --- a/Catalog/C6_fits_Catalog.py +++ b/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., diff --git a/Catalog/FGS_Catalog.py b/Catalog/FGS_Catalog.py index b1e57afc57a9561cc4d37f4a090f3460ca99246d..3b0c23024d87a7e2cc8074692a3a401657f78e55 100644 --- a/Catalog/FGS_Catalog.py +++ b/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., diff --git a/Catalog/NGPCatalog.py b/Catalog/NGPCatalog.py index 73433f72a1b7732033b0597c912fb4b147d41aa7..f37b73c467297ba5f514e2d374640b25b5626175 100644 --- a/Catalog/NGPCatalog.py +++ b/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) diff --git a/Catalog/fd_test_C6.py b/Catalog/fd_test_C6.py new file mode 100644 index 0000000000000000000000000000000000000000..0b7909a3e3c7dfe3dce5cc2c3d0f55e66a9011a9 --- /dev/null +++ b/Catalog/fd_test_C6.py @@ -0,0 +1,524 @@ +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 diff --git a/Catalog/testCat_galaxy.py b/Catalog/testCat_galaxy.py index 4859bad612fc4fe23120cc77e7aad3a533d90247..8ecd698e421cadd9aeb742701cab191d701776c8 100644 --- a/Catalog/testCat_galaxy.py +++ b/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., diff --git a/Catalog/wcs_test_C6.py b/Catalog/wcs_test_C6.py new file mode 100644 index 0000000000000000000000000000000000000000..896ab4ae5e4cd86a861509dc5556f709c7769e12 --- /dev/null +++ b/Catalog/wcs_test_C6.py @@ -0,0 +1,520 @@ +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 diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 6959253b38494409806d41217bb27e44e77144cb..8bf9d28c8c00e18f4c3d65ae43eecab5b35e469e 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/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] diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 57984092efe1af7aa98f34eaec3d02720f99fc9b..cc0b4ce81623d1b551952b951eda084c43a6aca7 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/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: @@ -552,6 +538,21 @@ class Chip(FocalPlane): else: 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: diff --git a/ObservationSim/Instrument/FilterParam.py b/ObservationSim/Instrument/FilterParam.py index 80057a52a157963efa260632cd77c48bd4a530e2..8f87b8fffd1f661b7c1eb2ef3eec352a5e24f94f 100755 --- a/ObservationSim/Instrument/FilterParam.py +++ b/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], diff --git a/ObservationSim/Instrument/FocalPlane.py b/ObservationSim/Instrument/FocalPlane.py index 55ef478911366aefdc433a6cae340206eb9ff614..ddca78f7aae07b6c24d914d3b36ba7dc58c37bd8 100755 --- a/ObservationSim/Instrument/FocalPlane.py +++ b/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) diff --git a/ObservationSim/Instrument/data/ccd/.DS_Store b/ObservationSim/Instrument/data/ccd/.DS_Store deleted file mode 100644 index f5986581377d6596ff3f73a930a273a694ebd847..0000000000000000000000000000000000000000 Binary files a/ObservationSim/Instrument/data/ccd/.DS_Store and /dev/null differ diff --git a/ObservationSim/Instrument/data/field_distortion/FieldDistModel_v2.0.pickle b/ObservationSim/Instrument/data/field_distortion/FieldDistModel_v2.0.pickle new file mode 100644 index 0000000000000000000000000000000000000000..7614524fd9e175295f4a4bf2001c257d5f2c6a84 Binary files /dev/null and b/ObservationSim/Instrument/data/field_distortion/FieldDistModel_v2.0.pickle differ diff --git a/ObservationSim/Instrument/data/field_distortion/__pycache__/__init__.cpython-39.pyc b/ObservationSim/Instrument/data/field_distortion/__pycache__/__init__.cpython-39.pyc deleted file mode 100644 index 347511c2b740a42e0dff664b5c77a40f93200323..0000000000000000000000000000000000000000 Binary files a/ObservationSim/Instrument/data/field_distortion/__pycache__/__init__.cpython-39.pyc and /dev/null differ diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py index a21edbaf4e8dd561ba28cd71cdbdfced98602eab..a5e7e6b84e2cb4a6970b6f0af01c4aad8884b386 100644 --- a/ObservationSim/MockObject/CatalogBase.py +++ b/ObservationSim/MockObject/CatalogBase.py @@ -1,6 +1,7 @@ import numpy as np import galsim import copy +import cmath from astropy.table import Table from abc import abstractmethod, ABCMeta @@ -72,6 +73,16 @@ class CatalogBase(metaclass=ABCMeta): "parallax":1e-9 } 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): diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 3db13dcdbec118746a97e918e488bd0b0dbbd3be..6bda70164cb33d8585a40aff447f1df1f471fb3f 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/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,7 +49,8 @@ class Galaxy(MockObject): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return -1 for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -59,7 +58,8 @@ class Galaxy(MockObject): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return -1 ratio = sub/full @@ -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,7 +98,8 @@ class Galaxy(MockObject): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return 2, None nphotons_sum = 0 @@ -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,7 +152,8 @@ class Galaxy(MockObject): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(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_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk + gal_temp = gal_temp.shear(gal_shear) - gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk + 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,7 +229,8 @@ class Galaxy(MockObject): else: # Return code 0: object photons missed this detector print("obj %s missed"%(self.id)) - self.logger.info("obj %s missed"%(self.id)) + if self.logger: + self.logger.info("obj %s missed"%(self.id)) return 0, pos_shear # # [C6 TEST] @@ -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 diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 0b5ff7fc269577841ac64a043b7481ef49001644..ae7b7081651d4aa77857f74e859bca7360b39d11 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/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): @@ -79,8 +80,10 @@ class MockObject(object): dx = x - self.x_nominal 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,7 +135,8 @@ class MockObject(object): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return 2, None nphotons_sum = 0 @@ -164,7 +168,8 @@ class MockObject(object): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) # return False continue @@ -214,7 +219,8 @@ class MockObject(object): else: # Return code 0: object photons missed this detector print("obj %s missed"%(self.id)) - self.logger.info("obj %s missed"%(self.id)) + if self.logger: + self.logger.info("obj %s missed"%(self.id)) return 0, pos_shear del photons_list diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so deleted file mode 100755 index df39d4dd7a3e9e91461e47547cdda61b0da5a58c..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so deleted file mode 100755 index 8f65db26166811166f208b031d4d38bcd51bc9e6..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.7/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so deleted file mode 100755 index 2ed89c7e7d08aa2449d0d8aad1f2add6cb954021..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/interp.cpython-38-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/interp.cpython-38-x86_64-linux-gnu.so deleted file mode 100755 index 355c81e7e32f056a1241cd702cc401c41ae88800..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.linux-x86_64-3.8/disperse_c/interp.cpython-38-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/disperse.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/disperse.cpython-37m-darwin.so deleted file mode 100755 index 5197cdcf8cae7a9ff2d15749ca310ef4a064f678..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/disperse.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/interp.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/interp.cpython-37m-darwin.so deleted file mode 100755 index d359bb3ccb84e0bb176179d1c2978df9be128f6d..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/MockObject/disperse_c/interp.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/disperse.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/disperse.cpython-37m-darwin.so deleted file mode 100755 index 4d1d315c5ca3555585365a486b413f892ef9f152..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/disperse.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/interp.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/interp.cpython-37m-darwin.so deleted file mode 100755 index 0061b7af9b20c163a478bf13b77b6ecb43ee0d2a..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/lib.macosx-10.9-x86_64-3.7/disperse_c/interp.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/disperse.o b/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/disperse.o deleted file mode 100644 index c99c2b26bfc359149d028fd8017ed40207b65f40..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/disperse.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/interp.o b/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/interp.o deleted file mode 100644 index 743b3f6ab802a6d5418696d6a3e677ad7002aea4..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.7/disperse_c/interp.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/disperse.o b/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/disperse.o deleted file mode 100644 index 544f093bb01d078415c9b07d2289895615a2b802..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/disperse.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/interp.o b/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/interp.o deleted file mode 100644 index 0524e54b0594ae059d56841ccc43cdd1ddbdafac..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.linux-x86_64-3.8/disperse_c/interp.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/disperse.o b/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/disperse.o deleted file mode 100644 index e72ab2785974393e518eb6af609db30cf4777668..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/disperse.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/interp.o b/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/interp.o deleted file mode 100644 index a2423056cb36c54b24cb2d1163597ded1c57adbb..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/build/temp.macosx-10.9-x86_64-3.7/disperse_c/interp.o and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/.DS_Store b/ObservationSim/MockObject/SpecDisperser/disperse_c/.DS_Store deleted file mode 100644 index 927866aaf638fb901d10c70bd308faa67c730ec7..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/.DS_Store and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-darwin.so deleted file mode 100755 index 4d1d315c5ca3555585365a486b413f892ef9f152..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so deleted file mode 100755 index df39d4dd7a3e9e91461e47547cdda61b0da5a58c..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-37m-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-darwin.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-darwin.so deleted file mode 100755 index 0061b7af9b20c163a478bf13b77b6ecb43ee0d2a..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-darwin.so and /dev/null differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so deleted file mode 100755 index 8f65db26166811166f208b031d4d38bcd51bc9e6..0000000000000000000000000000000000000000 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-37m-x86_64-linux-gnu.so and /dev/null differ diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 0822789e40b974826e0b836b79e187c6f5b80787..a70f99788a1f5a955072670dcbd17ce14defa4e4 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -48,7 +48,7 @@ class Observation(object): self.filter_list.append(filt) self.all_filter.append(filt) - def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, shear_cat_file=None, cat_dir=None, sed_dir=None): + def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, cat_dir=None, sed_dir=None): chip_output.Log_info(':::::::::::::::::::Current Pointing Information::::::::::::::::::') chip_output.Log_info("RA: %f, DEC; %f" % (pointing.ra, pointing.dec)) @@ -68,8 +68,7 @@ class Observation(object): chip_output.Log_error("unrecognized PSF model type!!", flush=True) # Figure out shear fields - if shear_cat_file is not None: - self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config, shear_cat_file=shear_cat_file) + self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config) # Apply astrometric simulation for pointing if self.config["obs_setting"]["enable_astrometric_model"]: @@ -111,8 +110,8 @@ class Observation(object): if self.config["obs_setting"]["enable_straylight_model"]: filt.setFilterStrayLightPixel(jtime = pointing.jdt, sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), pointing_radec = np.array([pointing.ra,pointing.dec]), sun_pos = np.array([pointing.sun_x,pointing.sun_y,pointing.sun_z])) - print("========================sky pix========================") - print(filt.sky_background) + chip_output.Log_info("========================sky pix========================") + chip_output.Log_info(filt.sky_background) if chip.survey_type == "photometric": sky_map = None @@ -160,7 +159,7 @@ class Observation(object): cut_filter = temp_filter if self.config["ins_effects"]["field_dist"] == True: - self.fd_model = FieldDistortion(chip=chip) + self.fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) else: self.fd_model = None @@ -188,6 +187,8 @@ class Observation(object): timestamp = pointing.timestamp, exptime = pointing.exp_time, readoutTime = 40.) + + chip_wcs = galsim.FitsWCS(header=h_ext) for j in range(self.nobj): @@ -243,12 +244,6 @@ class Observation(object): obj.g1, obj.g2 = 0., 0. else: obj.g1, obj.g2 = self.g1_field, self.g2_field - elif self.config["shear_setting"]["shear_type"] == "extra": - try: - # [TODO]: every object with individual shear from input catalog(s) - obj.g1, obj.g2 = self.g1_field[j], self.g2_field[j] - except: - chip_output.Log_error("failed to load external shear.") elif self.config["shear_setting"]["shear_type"] == "catalog": pass else: @@ -256,7 +251,7 @@ class Observation(object): raise ValueError("Unknown shear input") # Get position of object on the focal plane - pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=h_ext) + pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=h_ext) # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # Otherwise they will be considered missed objects @@ -410,7 +405,7 @@ class Observation(object): chip_output.Log_info("check running:2: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) - def runExposure_MPI_PointingList(self, pointing_list, shear_cat_file=None, chips=None, use_mpi=False): + def runExposure_MPI_PointingList(self, pointing_list,chips=None, use_mpi=False): if use_mpi: comm = MPI.COMM_WORLD ind_thread = comm.Get_rank() diff --git a/ObservationSim/PSF/FieldDistortion.py b/ObservationSim/PSF/FieldDistortion.py index 95dc40e7547e57e8e1545c42ac6444708271e09b..7c95fd6b7b3369cb2ae5792569dc293989d02a6c 100644 --- a/ObservationSim/PSF/FieldDistortion.py +++ b/ObservationSim/PSF/FieldDistortion.py @@ -1,9 +1,10 @@ import galsim import numpy as np +import cmath class FieldDistortion(object): - def __init__(self, chip, fdModel=None, fdModel_path=None): + def __init__(self, chip, fdModel=None, fdModel_path=None, img_rot=0.): if fdModel is None: if hasattr(chip, 'fdModel'): self.fdModel = chip.fdModel @@ -15,6 +16,7 @@ class FieldDistortion(object): raise ValueError("Error: no field distortion model has been specified!") else: self.fdModel = fdModel + self.img_rot = img_rot self.ifdModel = self.fdModel["wave1"] self.ixfdModel = self.ifdModel["xImagePos"] self.iyfdModel = self.ifdModel["yImagePos"] @@ -42,7 +44,7 @@ class FieldDistortion(object): return False return True - def get_distorted(self, chip, pos_img, bandpass=None): + def get_distorted(self, chip, pos_img, bandpass=None, img_rot=None): """ Get the distored position for an undistorted image position Parameters: @@ -58,14 +60,14 @@ class FieldDistortion(object): """ if not self.isContainObj_FD(chip=chip, pos_img=pos_img): return galsim.PositionD(-1, -1), None + if not img_rot: + img_rot = np.radians(self.img_rot) + else: + img_rot = np.radians(img_rot) x, y = pos_img.x, pos_img.y x = self.ixfdModel(x, y)[0][0] y = self.iyfdModel(x, y)[0][0] - ix_dx = self.ifx_dx(x, y) - ix_dy = self.ifx_dy(x, y) - iy_dx = self.ify_dx(x, y) - iy_dy = self.ify_dy(x, y) - if self.irsModel is not None: + if self.irsModel: # x1LowI, x1UpI, y1LowI, y1UpI = self.irsModel["interpLimit"] # if (x1LowI-x)*(x1UpI-x) <=0 and (y1LowI-y)*(y1UpI-y)<=0: # dx = self.ixrsModel(x, y)[0][0] @@ -88,8 +90,21 @@ class FieldDistortion(object): ix_dy = self.ifx_dy(x, y) + self.irx_dy(x, y) iy_dx = self.ify_dx(x, y) + self.iry_dx(x, y) iy_dy = self.ify_dy(x, y) + self.iry_dy(x, y) + else: + ix_dx = self.ifx_dx(x, y) + ix_dy = self.ifx_dy(x, y) + iy_dx = self.ify_dx(x, y) + iy_dy = self.ify_dy(x, y) g1k_fd = 0.0 + (iy_dy - ix_dx) / (iy_dy + ix_dx) g2k_fd = 0.0 - (iy_dx + ix_dy) / (iy_dy + ix_dx) + + # [TODO] [TESTING] Rotate the shear: + g_abs = np.sqrt(g1k_fd**2 + g2k_fd**2) + phi = cmath.phase(complex(g1k_fd, g2k_fd)) + # g_abs = 0.7 + g1k_fd = g_abs * np.cos(phi + 2*img_rot) + g2k_fd = g_abs * np.sin(phi + 2*img_rot) + fd_shear = galsim.Shear(g1=g1k_fd, g2=g2k_fd) return galsim.PositionD(x, y), fd_shear diff --git a/ObservationSim/_util.py b/ObservationSim/_util.py index 6c5af297aaab739f2f05d7dbd6899e11b8023a27..9f12f5d2c03feea2fba56cb2ad6e7207d08a7167 100755 --- a/ObservationSim/_util.py +++ b/ObservationSim/_util.py @@ -144,8 +144,8 @@ def makeSubDir_PointingList(path_dict, config, pointing_ID=0): pass return subImgdir, prefix -def get_shear_field(config, shear_cat_file=None): - if not config["shear_setting"]["shear_type"] in ["constant", "extra", "catalog"]: +def get_shear_field(config): + if not config["shear_setting"]["shear_type"] in ["constant", "catalog"]: raise ValueError("Please set a right 'shear_method' parameter.") if config["shear_setting"]["shear_type"] == "constant": @@ -153,18 +153,6 @@ def get_shear_field(config, shear_cat_file=None): g2 = config["shear_setting"]["reduced_g2"] nshear = 1 # TODO logging - elif config["shear_setting"]["shear_type"] == "extra": - # TODO logging - if not os.path.exists(shear_cat_file): - raise ValueError("Cannot find external shear catalog file.") - try: - shearCat = np.loadtxt(shear_cat_file) - nshear = shearCat.shape[0] - g1, g2 = shearCat[:, 0], shearCat[:, 1] - except: - print("Failed to read the shear catalog file.") - print("Setting to no shear.") - g1, g2 = 0., 0. else: g1, g2 = 0., 0. nshear = 0 diff --git a/README.md b/README.md index 133bf4e1aa9c2a1bcd9b93c3c258529705ea1f9d..b050654cef82d157d6b9f902d354e4a9afb578cb 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ * 导星探测器仿真 * 仪器参数更新，bug修复等 * 2022.08.23: 更新至最新的透过率曲线，在图像头文件中加入版本号。仿真软件版本号进行重新标记，当前版本为1.0.4 -* 2022.06.20: 版本v0.5上线，请于release_v0.5 branch中进行clone操作(注：目前已重新设置为v1.0，release_v1.0.0) +* 2022.06.20: 版本v0.5上线，release_v0.5 branch(注：目前已重新设置为v1.0，release_v1.0.0) * 2022.04.13: 修复了天测模块中的调用错误 * 2022.04.08: 修复了仅输出星表选项("out_cat_only")存在的问题 diff --git a/config/config_50sqdeg.yaml b/config/config_50sqdeg.yaml new file mode 100644 index 0000000000000000000000000000000000000000..c0f4e7ffc4889323cb35726e64871bb6082d6005 --- /dev/null +++ b/config/config_50sqdeg.yaml @@ -0,0 +1,221 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# Base diretories and naming setup +# Can add some of the command-line arguments here as well; +# OK to pass either way or both, as long as they are consistent +work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +data_dir: "/share/simudata/CSSOSDataProductsSims/data/" +run_name: "rotate_0" + +# Whether to use MPI +run_option: + use_mpi: NO + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "" + # star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + star_cat: "starcat/" + galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" + # AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "qsocat/qsosed/" + # AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: YES + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + seed_Av: 121212 # Seed for generating random intrinsic extinction + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscopic": simulate slitless spectroscopic chips only + # "FGS": simulate FGS chips only (31-42) + # "All": simulate full focal plane + survey_type: "Photometric" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/share/home/fangyuedong/50sqdeg_pointings/" + pointing_file: "pointing_50_combined.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [0] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [8] + + # Whether to enable astrometric modeling + enable_astrometric_model: True + + # Whether to enable straylight model + enable_straylight_model: True + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: ON # Whether to add field distortions + add_back: ON # Whether to add sky background + add_dark: ON # Whether to add dark noise + add_readout: ON # Whether to add read-out (Gaussian) noise + add_bias: ON # Whether to add bias-level to images + bias_16channel: ON # Whether to add different biases for 16 channels + gain_16channel: ON # Whether to make different gains for 16 channels + shutter_effect: ON # Whether to add shutter effect + flat_fielding: ON # Whether to add flat-fielding effect + prnu_effect: ON # Whether to add PRNU effect + non_linear: ON # Whether to add non-linearity + cosmic_ray: ON # Whether to add cosmic-ray + cray_differ: ON # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: ON # Whether to simulate CTE trails + saturbloom: ON # Whether to simulate Saturation & Blooming + add_badcolumns: ON # Whether to add bad columns + add_hotpixels: ON # Whether to add hot pixels + add_deadpixels: ON # Whether to add dead(dark) pixels + bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file diff --git a/config/config_C6.yaml b/config/config_C6.yaml index 087de4f806734afba835509cb784ee3fe5dc05ef..c130e59aee1b6bfa299f5ce6d2d5807bd86eff46 100644 --- a/config/config_C6.yaml +++ b/config/config_C6.yaml @@ -9,13 +9,13 @@ # Base diretories and naming setup # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent -work_dir: "/share/home/fangyuedong/fgs_sim/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "C6_spectroscopic_20230225" +run_name: "profile_C6" # Whether to use MPI run_option: - use_mpi: YES + use_mpi: NO # NOTE: "n_threads" paramters is currently not used in the backend # simulation codes. It should be implemented later in the web frontend # in order to config the number of threads to request from NAOC cluster @@ -100,13 +100,13 @@ obs_setting: # - give a list of indexes of chips: [ip_1, ip_2...] # - run all chips: null # Note: for all pointings - run_chips: [24] + run_chips: [8] # Whether to enable astrometric modeling enable_astrometric_model: True # Whether to enable straylight model - enable_straylight_model: False + enable_straylight_model: True # Cut by saturation magnitude in which band? cut_in_band: "z" @@ -144,17 +144,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "catalog" # For constant shear filed reduced_g1: 0. reduced_g2: 0. - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### @@ -162,25 +157,25 @@ ins_effects: # switches # Note: bias_16channel, gain_16channel, and shutter_effect # is currently not applicable to "FGS" observations - field_dist: ON # Whether to add field distortions - add_back: ON # Whether to add sky background - add_dark: ON # Whether to add dark noise - add_readout: ON # Whether to add read-out (Gaussian) noise - add_bias: ON # Whether to add bias-level to images - bias_16channel: ON # Whether to add different biases for 16 channels - gain_16channel: ON # Whether to make different gains for 16 channels - shutter_effect: ON # Whether to add shutter effect - flat_fielding: ON # Whether to add flat-fielding effect - prnu_effect: ON # Whether to add PRNU effect - non_linear: ON # Whether to add non-linearity - cosmic_ray: ON # Whether to add cosmic-ray - cray_differ: ON # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: ON # Whether to simulate CTE trails - saturbloom: ON # Whether to simulate Saturation & Blooming - add_badcolumns: ON # Whether to add bad columns - add_hotpixels: ON # Whether to add hot pixels - add_deadpixels: ON # Whether to add dead(dark) pixels - bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect + field_dist: YES # Whether to add field distortions + add_back: YES # Whether to add sky background + add_dark: YES # Whether to add dark noise + add_readout: YES # Whether to add read-out (Gaussian) noise + add_bias: YES # Whether to add bias-level to images + bias_16channel: YES # Whether to add different biases for 16 channels + gain_16channel: YES # Whether to make different gains for 16 channels + shutter_effect: YES # Whether to add shutter effect + flat_fielding: YES # Whether to add flat-fielding effect + prnu_effect: YES # Whether to add PRNU effect + non_linear: YES # Whether to add non-linearity + cosmic_ray: YES # Whether to add cosmic-ray + cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: YES # Whether to simulate CTE trails + saturbloom: YES # Whether to simulate Saturation & Blooming + add_badcolumns: YES # Whether to add bad columns + add_hotpixels: YES # Whether to add hot pixels + add_deadpixels: YES # Whether to add dead(dark) pixels + bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect # Values: # default values have been defined individually for each chip in: diff --git a/config/config_C6_fits.yaml b/config/config_C6_fits.yaml index 4837152ada2e5bd955ad36cb88655b8f713366ff..074a6279f52ace03ad3dc5958c7f908b30099e87 100644 --- a/config/config_C6_fits.yaml +++ b/config/config_C6_fits.yaml @@ -147,17 +147,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "catalog" # For constant shear filed reduced_g1: 0. reduced_g2: 0. - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### diff --git a/config/config_C6_test_wcs.yaml b/config/config_C6_test_wcs.yaml new file mode 100644 index 0000000000000000000000000000000000000000..9aff707ca4ccbf14956897c3ee6dd53c0728e209 --- /dev/null +++ b/config/config_C6_test_wcs.yaml @@ -0,0 +1,220 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# Base diretories and naming setup +# Can add some of the command-line arguments here as well; +# OK to pass either way or both, as long as they are consistent +work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +data_dir: "/share/simudata/CSSOSDataProductsSims/data/" +run_name: "fd_shear_g2m" + +# Whether to use MPI +run_option: + use_mpi: YES + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "Catalog_C6_20221212" + star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + galaxy_cat: "cat2CSSTSim_bundle/" + AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "quickspeclib_ross13.fits" + AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: YES + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + seed_Av: 121212 # Seed for generating random intrinsic extinction + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscopic": simulate slitless spectroscopic chips only + # "FGS": simulate FGS chips only (31-42) + # "All": simulate full focal plane + survey_type: "All" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" + pointing_file: "pointing_radec_246.5_40.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [0] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [9] + + # Whether to enable astrometric modeling + enable_astrometric_model: False + + # Whether to enable straylight model + enable_straylight_model: False + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: YES # Whether to add field distortions + add_back: YES # Whether to add sky background + add_dark: YES # Whether to add dark noise + add_readout: YES # Whether to add read-out (Gaussian) noise + add_bias: YES # Whether to add bias-level to images + bias_16channel: OFF # Whether to add different biases for 16 channels + gain_16channel: OFF # Whether to make different gains for 16 channels + shutter_effect: OFF # Whether to add shutter effect + flat_fielding: OFF # Whether to add flat-fielding effect + prnu_effect: OFF # Whether to add PRNU effect + non_linear: OFF # Whether to add non-linearity + cosmic_ray: OFF # Whether to add cosmic-ray + cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: OFF # Whether to simulate CTE trails + saturbloom: OFF # Whether to simulate Saturation & Blooming + add_badcolumns: OFF # Whether to add bad columns + add_hotpixels: OFF # Whether to add hot pixels + add_deadpixels: OFF # Whether to add dead(dark) pixels + bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file diff --git a/config/config_NGP.yaml b/config/config_NGP.yaml index 37aab5a891c23582967de69908af5a46251b935f..2a986b37c2d949df41358c0fcc805668497f818f 100644 --- a/config/config_NGP.yaml +++ b/config/config_NGP.yaml @@ -138,17 +138,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "constant" # For constant shear filed reduced_g1: 0.026 reduced_g2: 0.015 - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### diff --git a/config/config_example.yaml b/config/config_example.yaml index a3cc4026bb662776dfa4defe8abf06e6175d21c1..5dcec7bdf13f0729233726deade29cef8bc86cd1 100644 --- a/config/config_example.yaml +++ b/config/config_example.yaml @@ -131,17 +131,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "constant" # For constant shear filed reduced_g1: 0.026 reduced_g2: 0.015 - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### diff --git a/config/config_fgs.yaml b/config/config_fgs.yaml index 00bca3b67c2b5e77526b3bca3928d064fe4eb09a..3a0d021da04ac8428696fe10ae2396e90bbe0c20 100644 --- a/config/config_fgs.yaml +++ b/config/config_fgs.yaml @@ -143,17 +143,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "catalog" # For constant shear filed reduced_g1: 0. reduced_g2: 0. - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### diff --git a/config/config_testCASE.yaml b/config/config_testCASE.yaml index 062853b86f506f120a092606d9e5cd4a4a673978..503d753a15be2031b5beb933fcd5e33c2b1c0af5 100644 --- a/config/config_testCASE.yaml +++ b/config/config_testCASE.yaml @@ -146,17 +146,12 @@ shear_setting: # Options to generate mock shear field: # "constant": all galaxies are assigned a constant reduced shear # "catalog": from catalog - # "extra": from seprate file shear_type: "catalog" # For constant shear filed reduced_g1: 0. reduced_g2: 0. - # Extra shear catalog - # (currently not used) - # shear_cat: "mockShear.cat" - ############################################### # Instrumental effects setting ############################################### diff --git a/config/test_fd_C6.yaml b/config/test_fd_C6.yaml new file mode 100644 index 0000000000000000000000000000000000000000..3ecc3feb6c06d6e5f36614d8855a40568939f8e1 --- /dev/null +++ b/config/test_fd_C6.yaml @@ -0,0 +1,220 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# Base diretories and naming setup +# Can add some of the command-line arguments here as well; +# OK to pass either way or both, as long as they are consistent +work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +data_dir: "/share/simudata/CSSOSDataProductsSims/data/" +run_name: "fd_shear_g2m" + +# Whether to use MPI +run_option: + use_mpi: NO + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "Catalog_C6_20221212" + star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + galaxy_cat: "cat2CSSTSim_bundle/" + AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "quickspeclib_ross13.fits" + AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: YES + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + seed_Av: 121212 # Seed for generating random intrinsic extinction + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscopic": simulate slitless spectroscopic chips only + # "FGS": simulate FGS chips only (31-42) + # "All": simulate full focal plane + survey_type: "Photometric" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" + pointing_file: "pointing_radec_246.5_40.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [0] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [9] + + # Whether to enable astrometric modeling + enable_astrometric_model: True + + # Whether to enable straylight model + enable_straylight_model: False + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: YES # Whether to add field distortions + add_back: YES # Whether to add sky background + add_dark: YES # Whether to add dark noise + add_readout: NO # Whether to add read-out (Gaussian) noise + add_bias: YES # Whether to add bias-level to images + bias_16channel: NO # Whether to add different biases for 16 channels + gain_16channel: NO # Whether to make different gains for 16 channels + shutter_effect: NO # Whether to add shutter effect + flat_fielding: YES # Whether to add flat-fielding effect + prnu_effect: NO # Whether to add PRNU effect + non_linear: NO # Whether to add non-linearity + cosmic_ray: NO # Whether to add cosmic-ray + cray_differ: NO # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: NO # Whether to simulate CTE trails + saturbloom: NO # Whether to simulate Saturation & Blooming + add_badcolumns: NO # Whether to add bad columns + add_hotpixels: NO # Whether to add hot pixels + add_deadpixels: NO # Whether to add dead(dark) pixels + bright_fatter: NO # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file diff --git a/profile_C6.sh b/profile_C6.sh new file mode 100755 index 0000000000000000000000000000000000000000..49c79694a11d0082b595dfecf68e03b5ccca03a2 --- /dev/null +++ b/profile_C6.sh @@ -0,0 +1,17 @@ +#!/bin/bash + +date + +python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/csst-simulation/run_sim.py \ + --config_file config_50sqdeg.yaml \ + --catalog C6_50sqdeg \ + -c /share/home/fangyuedong/csst-simulation/config + + # --config_file config_C6.yaml \ + # --catalog C6_Catalog \ + # -c /share/home/fangyuedong/csst-simulation/config + + # --config_file test_fd_C6.yaml \ + # --catalog fd_test_C6 \ + # --config_file config_C6_test_wcs.yaml \ + # --catalog wcs_test_C6 \ diff --git a/run_C6.pbs b/run_C6.pbs index 9fd867219deb03f388f44b297b5bd130f5b8d5c5..4d4da6e734b114519e93df4147d881a555473325 100755 --- a/run_C6.pbs +++ b/run_C6.pbs @@ -1,18 +1,23 @@ #!/bin/bash - -#PBS -N SIMS -#PBS -l nodes=wcl-1:ppn=60 -###PBS -l nodes=wcl-1:ppn=24+wcl-2:ppn=24+wcl-3:ppn=24+wcl-4:ppn=24+wcl-5:ppn=24+wcl-6:ppn=24 +#PBS -N C6_TEST +#PBS -l walltime=70:00:00 +#PBS -j oe +##PBS -l nodes=wcl-1:ppn=16+wcl-2:ppn=16 +#PBS -q batch #PBS -u fangyuedong -###PBS -j oe cd $PBS_O_WORKDIR -NP=60 +NP=96 +hostfile=wcl-1,wcl-2 date -mpirun -np $NP python3 /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ - --config_file config_C6.yaml \ - --catalog C6_Catalog \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config +mpirun --oversubscribe -H $hostfile -np $NP python /share/home/fangyuedong/csst-simulation/run_sim.py \ + --config_file config_50sqdeg.yaml \ + --catalog C6_50sqdeg \ + -c /share/home/fangyuedong/csst-simulation/config + # --config_file config_C6.yaml \ + # --catalog C6_Catalog \ + # -c /share/home/fangyuedong/csst-simulation/config + diff --git a/test_C6.sh b/test_C6.sh index 82b04e595b71553c1c138f9a1eb5194347677f0e..d976f42dba8cd3185ffeaa61062869a1fc90ee17 100755 --- a/test_C6.sh +++ b/test_C6.sh @@ -2,8 +2,8 @@ date -python3 /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ +python3 /share/home/fangyuedong/sim_v2/csst-simulation/run_sim.py \ --config_file config_C6.yaml \ --catalog C6_Catalog \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config + -c /share/home/fangyuedong/sim_v2/csst-simulation/config