diff --git a/catalog/C10_Catalog.py b/catalog/C10_Catalog.py new file mode 100644 index 0000000000000000000000000000000000000000..8bbebcb7ebcde5dea8677f4c143d29f92a431608 --- /dev/null +++ b/catalog/C10_Catalog.py @@ -0,0 +1,574 @@ +import os +import galsim +import random +import copy +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 observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, ExtinctionMW +from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist +from observation_sim.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 +import astropy.constants as atcons + +import ctypes + +try: + import importlib.resources as pkg_resources +except ImportError: + # Try backported to PY<37 'importlib_resources' + import importlib_resources as pkg_resources + +NSIDE = 128 + +pointing_file_list = ['pointing_C10_1000deg_lat80.dat', 'pointing_C10_1000deg_latm80_.dat', 'pointing_C10_1000deg_p180_m30.dat', 'pointing_C10_1000deg_p320_m40.dat'] +star_file_list = ['CSST_C10_1000sqd_1.hdf5', 'CSST_C10_1000sqd_2_Fornax.hdf5', 'CSST_C10_1000sqd_3.hdf5', 'CSST_C10_1000sqd_4.hdf5'] + + +def get_galaxy_qso_list(config): + cat_dir = config["catalog_options"]["input_path"]["cat_dir"] + with open(cat_dir+"galcat_C10/qsocat/gal_C10_file", 'r', encoding='utf-8') as fn1: + fn1_list = fn1.readlines() + bundle_file_list = [line.strip() for line in fn1_list] + with open(cat_dir+"galcat_C10/qsocat/qso_sed_C10_file", 'r', encoding='utf-8') as fn2: + fn2_list = fn2.readlines() + qsosed_file_list = [line.strip() for line in fn2_list] + return bundle_file_list, qsosed_file_list + + +class StarParm(ctypes.Structure): + _fields_ = [ + ('logte', ctypes.c_float), + ('logg', ctypes.c_float), + ('Mass', ctypes.c_float), + ('Av', ctypes.c_float), + ('mu0', ctypes.c_float), + ('Z', ctypes.c_float)] + + +def get_star_cat(config): + idx = pointing_file_list.index(os.path.basename(config['obs_setting']['pointing_file'])) + return_star_path = star_file_list[idx] + return return_star_path + + +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, config): + bundle_file_list, qsosed_file_list = get_galaxy_qso_list(config) + return qsosed_file_list[bundle_file_list.index(bundle_file_name)] + + +class Catalog(CatalogBase): + def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): + super().__init__() + self.cat_dir = config["catalog_options"]["input_path"]["cat_dir"] + + 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. + + # [TODO] Milky Way extinction + if "enable_mw_ext_gal" in config["catalog_options"] and config["catalog_options"]["enable_mw_ext_gal"]: + if "planck_ebv_map" not in config["catalog_options"]: + raise ValueError( + "Planck dust map must be given to enable Milky Way extinction calculation for galaxies.") + self.mw_ext = ExtinctionMW() + self.mw_ext.init_ext_model(model_name="odonnell") + self.mw_ext.load_Planck_ext( + file_path=config["catalog_options"]["planck_ebv_map"]) + + 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(self.config) + star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name) + self.star_path = os.path.join(self.cat_dir, star_path) + self.star_SED_path = config["catalog_options"]["SED_templates_path"]["star_SED"] + 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 = 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 = 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 = " av stellarmass dm teff logg feh" + self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " + + self.add_fmt = "%8.4f %8.4f %8.4f %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_output_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_star(self): + # self.tempSED_star = h5.File(self.star_SED_path,'r') + with pkg_resources.path('catalog.data', 'starSpecInterp.so') as ddl_path: + self.starDDL = ctypes.CDLL(str(ddl_path)) + self.starDDL.loadSpecLibs.argtypes = [ctypes.c_char_p, ctypes.c_char_p] + self.starDDL.loadExts.argtypes = [ctypes.c_char_p] + nwv = self.starDDL.loadSpecLibs(str.encode(os.path.join( + self.star_SED_path, 'file_BT-Settl_CSST_wl1000-24000_R1000.par')), str.encode(self.star_SED_path)) + self.starDDL.loadExts(str.encode(os.path.join( + self.star_SED_path, "Ext_odonnell94_R3.1_CSST_wl1000-24000_R1000.fits"))) + self.star_spec_len = nwv + + def _interp_star_sed(self, obj): + spec = (ctypes.c_float*self.star_spec_len)() + wave = (ctypes.c_float*self.star_spec_len)() + self.starDDL.interpSingleStar.argtypes = [ + ctypes.Structure, ctypes.POINTER(ctypes.c_float)] + # s=Star(obj.param['teff'], obj.param['grav''], obj.paramstar['mwmsc_mass'], obj.param['AV'], obj.param['DM'], obj.param['z_met']) + s = StarParm(obj.param['teff'], obj.param['logg'], obj.param['stellarMass'], + obj.param['av'], obj.param['DM'], obj.param['feh']) + + self.starDDL.interpSingleStar(s, spec, wave) + + rv_c = obj.param['rv']/(atcons.c.value/1000.) + Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c)) + wave_RV = wave*Doppler_factor + return wave_RV, np.power(10, spec[:]) + + 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 + ) + + # [TODO] get Milky Way extinction AVs + if "enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]: + MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + else: + MW_Av_arr = np.zeros(len(ra_arr)) + + 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] + + # [TODO] get Milky Way extinction AVs + param['mw_Av'] = MW_Av_arr[igals] + + if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): + continue + + # 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 + + # TEMP + self.ids += 1 + param['id'] = '%06d' % (int(pix_id)) + \ + '%06d' % (cat_id) + '%08d' % (igals) + + # Is this an Quasar? + param['qsoindex'] = gals['qsoindex'][igals] + if param['qsoindex'] == -1: + param['star'] = 0 # Galaxy + param['agnsed_file'] = "" + obj = Galaxy(param, logger=self.logger) + else: + param_qso = copy.deepcopy(param) + param_qso['star'] = 2 # Quasar + param_qso['agnsed_file'] = agnsed_file + # First add QSO model + obj = Quasar(param_qso, logger=self.logger) + # Need to deal with additional output columns + obj.additional_output_str = self.add_fmt % (0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0, 0.) + self.objs.append(obj) + # Then add host galaxy model + param['star'] = 0 # Galaxy + param['agnsed_file'] = "" + obj = Galaxy(param, logger=self.logger) + + # Need to deal with additional output columns for (host) galaxy + obj.additional_output_str = self.add_fmt % (0., 0., 0., 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['RA']) + # 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'] = '%06d' % (int(pix_id)) + '%08d' % (istars) + # param['sed_type'] = istars + # param['model_tag'] = '' + param['teff'] = stars['teff'][istars] + param['logg'] = stars['grav'][istars] + param['feh'] = stars['z_met'][istars] + param['stellarMass'] = stars['mass'][istars] + + param['av'] = stars['AV'][istars] + param['DM'] = stars['DM'][istars] + # param['z_met'] = stars['z_met'][istars] + + param['z'] = 0.0 + param['star'] = 1 # Star + + try: + obj = Star(param, logger=self.logger) + except Exception as e: + print(e) + + # Append additional output columns to the .cat file + obj.additional_output_str = self.add_fmt % (param["av"], param['stellarMass'], param['DM'], 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')['star_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_C10_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, self.config) + 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'] + # ) + wave, flux = self._interp_star_sed(obj) + 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 * (1.0 + obj.z) + else: + raise ValueError("Object type not known") + speci = interpolate.interp1d(wave, flux, fill_value=0., bounds_error=False) + lamb = np.arange(2000, 11001+0.5, 0.5) + y = speci(lamb) + + # [TODO] Apply Milky Way extinction + if obj.type != 'star' and ("enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]): + y = self.mw_ext.apply_extinction(y, Av=obj.mw_Av) + + # 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'])) + + # integrate to get the magnitudes + if obj.type == 'quasar' or obj.type == 'galaxy': + 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/data/starSpecInterp.so b/catalog/data/starSpecInterp.so old mode 100755 new mode 100644 index 150d42c5d743164a22eef536f846317e85aae387..11899b46080fb4434dc6b8a20c89a23fd80d94dd Binary files a/catalog/data/starSpecInterp.so and b/catalog/data/starSpecInterp.so differ diff --git a/config/obs_config_SCI.yaml b/config/obs_config_SCI.yaml index fecdf01a758f2063adc778604daca624b7cb5287..fb8f538b2bbdfd3b42ec0cacb5e6c3538a014fa0 100644 --- a/config/obs_config_SCI.yaml +++ b/config/obs_config_SCI.yaml @@ -62,10 +62,10 @@ call_sequence: hot_pixels: YES dead_pixels: YES bad_columns: YES - # Apply response nonlinearity - nonlinearity: {} # Apply CCD Saturation & Blooming blooming: {} + # Apply response nonlinearity + nonlinearity: {} # Run CTE simulation # CTE_effect: {} # Add prescan and overscan diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index 523406e56d2161459015febfa7abe3b4e1e70d26..dbad1fbe44d4597d67d4f3532705d40d2704753e 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -154,13 +154,12 @@ class Observation(object): chip_output.Log_error(e) chip_output.Log_error("Failed simulation on step: %s" % (step)) break - - chip_output.Log_info("check running:1: 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))) del chip.img del chip.flat_img del chip.prnu_img del chip.shutter_img + chip_output.Log_info("check running:1: 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, chips=None): comm = MPI.COMM_WORLD @@ -239,5 +238,6 @@ class Observation(object): "finished running chip#%d..." % (chip.chipID)) for handler in chip_output.logger.handlers[:]: chip_output.logger.removeHandler(handler) + del chip_output gc.collect() process_counter += nchips_per_fp diff --git a/observation_sim/instruments/chip/libBF/diffusion_X1.c b/observation_sim/instruments/chip/libBF/diffusion_X1.c index 199b48cc94540419cc206dbf327b40aa13b2f849..b0259e8169130873d1bd3a2ec3fdcd08cf1fcdda 100644 --- a/observation_sim/instruments/chip/libBF/diffusion_X1.c +++ b/observation_sim/instruments/chip/libBF/diffusion_X1.c @@ -37,7 +37,7 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi { printf("Adding BF effect...\n"); //setup BF correlation fliter - float neX; + float neX, neXtemp; float neP1 = 50000; float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302}; float neP2 = 10000; @@ -56,15 +56,18 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi neX = arr_ima[j+i*ny]; if(neX >= 10000) { + neXtemp = neX; + if(neXtemp > 100000) + neXtemp = 100000; bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0; - bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575; - bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652; - bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335; - bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]); - bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118; - bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]); - bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083; - bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043; + bfa[0][1]=bfa[0][-1]=linearInterp(neXtemp, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575; + bfa[-1][0]=bfa[1][0]=linearInterp(neXtemp, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652; + bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neXtemp, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335; + bfa[0][-2]=bfa[0][2]=linearInterp(neXtemp, neP1, bfaP1[4], neP2, bfaP2[4]); + bfa[-2][0]=bfa[2][0]=linearInterp(neXtemp, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118; + bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neXtemp, neP1, bfaP1[6], neP2, bfaP2[6]); + bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neXtemp, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083; + bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neXtemp, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043; } else { diff --git a/observation_sim/mock_objects/Galaxy.py b/observation_sim/mock_objects/Galaxy.py index 03df3ebd84c35977f2421a8dd7dcdb4d5981ab5b..6575035695a25fe2bdc5261b259c8ec6c63104c5 100755 --- a/observation_sim/mock_objects/Galaxy.py +++ b/observation_sim/mock_objects/Galaxy.py @@ -115,7 +115,7 @@ class Galaxy(MockObject): # Set Galsim Parameters if self.getMagFilter(filt) <= 15 and (not big_galaxy): - folding_threshold = 5.e-4 + folding_threshold = 5.e-8 else: folding_threshold = 5.e-3 gsp = galsim.GSParams(folding_threshold=folding_threshold) @@ -170,8 +170,11 @@ class Galaxy(MockObject): # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) # Get PSF model + EXTRA = False + if self.getMagFilter(filt) <= filt.mag_saturation-1.: + EXTRA = True psf, pos_shear = psf_model.get_PSF( - chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) + chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold, extrapolate=EXTRA) if self.bfrac == 0: gal_temp = disk @@ -198,7 +201,7 @@ class Galaxy(MockObject): # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots # stamp = gal.drawImage(wcs=chip_wcs_local, method='phot', offset=offset, save_photons=True) - stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset) + stamp = gal.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset) if np.sum(np.isnan(stamp.array)) > 0: # ERROR happens return 2, pos_shear diff --git a/observation_sim/mock_objects/MockObject.py b/observation_sim/mock_objects/MockObject.py index 38608fa6ad8062350c89305ea600f241808e7153..8ef57be387eaeaf800b95bdc56a4de7ab2dec10d 100755 --- a/observation_sim/mock_objects/MockObject.py +++ b/observation_sim/mock_objects/MockObject.py @@ -122,7 +122,7 @@ class MockObject(object): return 2, None # Set Galsim Parameters if self.getMagFilter(filt) <= 15: - folding_threshold = 5.e-4 + folding_threshold = 5.e-8 else: folding_threshold = 5.e-3 gsp = galsim.GSParams(folding_threshold=folding_threshold) @@ -160,14 +160,17 @@ class MockObject(object): # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) # Get PSF model + EXTRA = False + if self.getMagFilter(filt) <= filt.mag_saturation-1.: + EXTRA = True psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, - folding_threshold=folding_threshold) + folding_threshold=folding_threshold, extrapolate=EXTRA) # star = galsim.DeltaFunction(gsparams=gsp) # star = star.withFlux(nphotons) # star = galsim.Convolve(psf, star) star = psf.withFlux(nphotons) - stamp = star.drawImage(wcs=chip_wcs_local, offset=offset) + stamp = star.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset) if np.sum(np.isnan(stamp.array)) > 0: continue stamp.setCenter(x_nominal, y_nominal) diff --git a/observation_sim/mock_objects/Stamp.py b/observation_sim/mock_objects/Stamp.py index d39ace7ec70edac6904dc2e992bb067a3c6c3872..c16c26f397d57f2548ef80c7edca32125eba5997 100644 --- a/observation_sim/mock_objects/Stamp.py +++ b/observation_sim/mock_objects/Stamp.py @@ -92,7 +92,7 @@ class Stamp(MockObject): else: gal = gal + gal_temp - stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset) + stamp = gal.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset) if np.sum(np.isnan(stamp.array)) > 0: # ERROR happens return 2, pos_shear diff --git a/observation_sim/psf/PSFInterp.py b/observation_sim/psf/PSFInterp.py index ab902d96859bdbcbe4612871ceac16b82b582abe..d21e75fcc81d030f41022da26a40e8e460ddc809 100644 --- a/observation_sim/psf/PSFInterp.py +++ b/observation_sim/psf/PSFInterp.py @@ -13,6 +13,7 @@ import galsim import h5py from observation_sim.psf.PSFModel import PSFModel +from observation_sim.psf._util import psf_extrapolate NPSF = 900 # ***# 30*30 @@ -324,7 +325,7 @@ class PSFInterp(PSFModel): return twave return -1 - def get_PSF(self, chip, pos_img, bandpass, galsimGSObject=True, findNeighMode='treeFind', folding_threshold=5.e-3, pointing_pa=0.0): + def get_PSF(self, chip, pos_img, bandpass, galsimGSObject=True, findNeighMode='treeFind', folding_threshold=5.e-3, pointing_pa=0.0, extrapolate=False, ngg=2048): """ Get the PSF at a given image position @@ -358,20 +359,18 @@ class PSFInterp(PSFModel): imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True) - ''' - ############TEST: START - TestGaussian = False - if TestGaussian: - gsx = galsim.Gaussian(sigma=0.04) - #pointing_pa = -23.433333 - imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array - ############TEST: END - ''' + if extrapolate is True: + ccdList = [6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23, 24, 25] + rr_trim_list = [72, 64, 96, 88, 64, 72, 72, 76, 72, 72, 76, 72, 72, 64, 88, 96, 64, 72] + imPSF = psf_extrapolate(imPSF, rr_trim=rr_trim_list[ccdList.index(chip.chipID)], ngg=ngg) if galsimGSObject: - imPSFt = np.zeros([257, 257]) - imPSFt[0:256, 0:256] = imPSF - # imPSFt[120:130, 0:256] = 1. + if extrapolate is True: + imPSFt = np.zeros([ngg+1, ngg+1]) + imPSFt[:-1, :-1] = imPSF + else: + imPSFt = np.zeros([257, 257]) + imPSFt[0:256, 0:256] = imPSF img = galsim.ImageF(imPSFt, scale=pixSize) gsp = galsim.GSParams(folding_threshold=folding_threshold) diff --git a/observation_sim/psf/_util.py b/observation_sim/psf/_util.py new file mode 100644 index 0000000000000000000000000000000000000000..45e1c02fcbf9ce3da871d5519eb7ee07f4e2d9ba --- /dev/null +++ b/observation_sim/psf/_util.py @@ -0,0 +1,65 @@ +import numpy as np +from scipy.interpolate import interp1d + + +def binningPSF(img, ngg): + imgX = img.reshape(ngg, img.shape[0]//ngg, ngg, img.shape[1]//ngg).mean(-1).mean(1) + return imgX + + +def radial_average_at_pixel(image, center_x, center_y, dr=10): + # Get coordinates relative to the specified center pixel (x, y) + y, x = np.indices(image.shape) + r = np.sqrt((x - center_x)**2 + (y - center_y)**2) + + # Set up bins + max_radius = int(r.max()) # Maximum distance from the center pixel + radial_bins = np.arange(0, max_radius, dr) + + # Compute average value in each bin + radial_means = [] + for i in range(len(radial_bins) - 1): + mask = (r >= radial_bins[i]) & (r < radial_bins[i + 1]) + if np.any(mask): + radial_means.append(image[mask].mean()) + else: + radial_means.append(0) # In case no pixels are in the bin + return radial_bins[:-1], radial_means # Exclude last bin since no mean calculated + + +def psf_extrapolate(psf, rr_trim=64, ngg=256): + # ngg = 256 + # extrapolate PSF + if True: + xim = np.arange(256)-128 + xim, yim = np.meshgrid(xim, xim) + rim = np.sqrt(xim**2 + yim**2) + + # rr_trim = 96 + psf_temp = psf + psf_temp[rim > rr_trim] = 0 + radii, means = radial_average_at_pixel(psf_temp, 128, 128, dr=4) + + radii_log = np.log(radii[1:]) + means_log = np.log(means[1:]) + + finite_mask = np.isfinite(means_log) + f_interp = interp1d(radii_log[finite_mask][:-1], means_log[finite_mask][:-1], kind='linear', fill_value="extrapolate") + + # ngg = 1024 + xim = np.arange(ngg)-int(ngg/2) + xim, yim = np.meshgrid(xim, xim) + rim = np.sqrt(xim**2 + yim**2) + rim[int(ngg/2), int(ngg/2)] = np.finfo(float).eps # 1e-7 + rim_log = np.log(rim) + y_new = f_interp(rim_log) + + arr = np.zeros([ngg, ngg]) + arr[int(ngg/2-128):int(ngg/2+128), int(ngg/2-128):int(ngg/2+128)] = np.log(psf_temp + np.finfo(float).eps) + arr[rim > rr_trim] = 0 + arr[arr == 0] = y_new[arr == 0] + psf = np.exp(arr) + psf[rim > int(ngg/2)] = 0 + imPSF = psf # binningPSF(psf, int(ngg/2)) + imPSF = imPSF/np.nansum(imPSF) + return imPSF diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 9e6177c416e1583679d7a4bad3648d3719b7accd..81cb5892128f0be3870d57089b2949b1de4485a1 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -238,11 +238,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): self.chip_output.Log_error("obj id: %s" % (obj.param['id'])) self.chip_output.Log_error(" e1: %.5f\n e2: %.5f\n size: %f\n bfrac: %f\n detA: %f\n g1: %.5f\n g2: %.5f\n" % ( obj.param['e1'], obj.param['e2'], obj.param['size'], obj.param['bfrac'], obj.param['detA'], obj.param['g1'], obj.param['g2'])) - # Unload SED: obj.unload_SED() del obj # gc.collect() + cat.starDDL.freeGlobeData() + del cat.starDDL if chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"] and chip.slsPSFOptim: # from observation_sim.instruments.chip import chip_utils as chip_utils diff --git a/tools/setConfig/reset_obs.py b/tools/setConfig/reset_obs.py new file mode 100644 index 0000000000000000000000000000000000000000..f05b9d2448458b35e6eb5c280d6eece1384d7991 --- /dev/null +++ b/tools/setConfig/reset_obs.py @@ -0,0 +1,122 @@ +from flask import Flask, render_template, request, redirect +import yaml +import ast + +app = Flask(__name__) + + +key_type_map = { + 'obs_type': str, + 'obs_type_code': str, + 'obs_id': str, + 'run_chips': list, + 'call_sequence': { + 'scie_obs': { + 'shutter_effect': bool, + 'flat_fielding': bool, + 'field_dist': bool, + }, + 'sky_background': { + 'shutter_effect': bool, + 'flat_fielding': bool, + 'enable_straylight_model': bool, + 'flat_level': None, + 'flat_level_filt': None, + }, + 'PRNU_effect': {}, + 'cosmic_rays': { + 'save_cosmic_img': bool, + }, + 'poisson_and_dark': { + 'add_dark': bool, + }, + 'bright_fatter': {}, + 'detector_defects': { + 'hot_pixels': bool, + 'dead_pixels': bool, + 'bad_columns': bool, + }, + 'nonlinearity': {}, + 'blooming': {}, + 'prescan_overscan': { + 'add_dark': bool, + }, + 'bias': { + 'bias_16channel': bool, + }, + 'readout_noise': {}, + 'gain': { + 'gain_16channel': bool, + }, + 'quantization_and_output': { + 'format_output': bool, + }, + }, +} + + +def convert_dict_values(d, key_type_map): + for key, value in d.items(): + if isinstance(value, dict): + convert_dict_values(value, key_type_map[key]) + elif key in key_type_map: + if key_type_map[key] is int: + d[key] = int(value) + if key_type_map[key] is float: + d[key] = float(value) + if key_type_map[key] is bool: + if d[key].lower() == 'yes' or d[key].lower() == 'true': + d[key] = True + else: + d[key] = False + if key_type_map[key] is str: + d[key] = str(value) + if key_type_map[key] is list: + d[key] = ast.literal_eval(value) + if key_type_map[key] is None: + d[key] = None + + +def load_yaml(): + with open('templates/obs_config_SCI.yaml', 'r') as file: + return yaml.safe_load(file) + + +def save_yaml(data): + convert_dict_values(data, key_type_map) + with open('config_reset/obs_config_SCI_reset.yaml', 'w') as file: + yaml.dump(data, file, default_flow_style=False, sort_keys=False) + + +def render_form(data, parent_key=''): + form_html = '' + for key, value in data.items(): + full_key = f"{parent_key}.{key}" if parent_key else key + if isinstance(value, dict): # 处理字典 + form_html += f"

{key}

{render_form(value, full_key)}
" + else: + form_html += f"" + form_html += f"
" + return form_html + + +@app.route('/', methods=['GET', 'POST']) +def index(): + if request.method == 'POST': + data = load_yaml() + for key in request.form: + keys = key.split('.') + temp = data + for k in keys[:-1]: + temp = temp[k] + temp[keys[-1]] = request.form[key] + save_yaml(data) + return redirect('/') + + data = load_yaml() + form_html = render_form(data) + return render_template('index_obs.html', form_html=form_html) + + +if __name__ == '__main__': + app.run(debug=True) diff --git a/tools/setConfig/reset_overall.py b/tools/setConfig/reset_overall.py new file mode 100644 index 0000000000000000000000000000000000000000..cacb2049210ed85717c103baa26fc2dc01be32d2 --- /dev/null +++ b/tools/setConfig/reset_overall.py @@ -0,0 +1,133 @@ +from flask import Flask, render_template, request, redirect +import yaml +import ast + +app = Flask(__name__) + +key_type_map = { + 'work_dir': str, + 'run_name': str, + 'project_cycle': int, + 'run_counter': int, + 'run_option': { + 'out_cat_only': bool, + }, + 'catalog_options': { + 'input_path': { + 'cat_dir': str, + 'star_cat': str, + 'galaxy_cat': str, + }, + 'SED_templates_path': { + 'star_SED': str, + 'galaxy_SED': str, + 'AGN_SED': str, + }, + 'star_only': bool, + 'galaxy_only': bool, + 'rotateEll': float, + 'enable_mw_ext_gal': bool, + 'planck_ebv_map': str, + }, + 'obs_setting': { + 'pointing_file': str, + 'obs_config_file': str, + 'run_pointings': list, + 'enable_astrometric_model': bool, + 'cut_in_band': str, + 'mag_sat_margin': float, + 'mag_lim_margin': float, + }, + 'psf_setting': { + 'psf_model': str, + 'psf_pho_dir': str, + 'psf_sls_dir': str, + }, + 'shear_setting': { + 'shear_type': str, + 'reduced_g1': float, + 'reduced_g2': float, + }, + 'output_setting': { + 'output_format': str, + 'shutter_output': bool, + 'prnu_output': bool, + }, + 'random_seeds': { + 'seed_poisson': int, + 'seed_CR': int, + 'seed_flat': int, + 'seed_prnu': int, + 'seed_gainNonUniform': int, + 'seed_biasNonUniform': int, + 'seed_rnNonUniform': int, + 'seed_badcolumns': int, + 'seed_defective': int, + 'seed_readout': int, + }, +} + + +def convert_dict_values(d, key_type_map): + for key, value in d.items(): + if isinstance(value, dict): + convert_dict_values(value, key_type_map[key]) + elif key in key_type_map: + if key_type_map[key] is int: + d[key] = int(value) + if key_type_map[key] is float: + d[key] = float(value) + if key_type_map[key] is bool: + if d[key].lower() == 'yes' or d[key].lower() == 'true': + d[key] = True + else: + d[key] = False + if key_type_map[key] is str: + d[key] = str(value) + if key_type_map[key] is list: + d[key] = ast.literal_eval(value) + + +def load_yaml(): + with open('templates/config_overall.yaml', 'r') as file: + return yaml.safe_load(file) + + +def save_yaml(data): + convert_dict_values(data, key_type_map) + with open('config_reset/config_overall_reset.yaml', 'w') as file: + yaml.dump(data, file, default_flow_style=False, sort_keys=False) + + +def render_form(data, parent_key=''): + form_html = '' + for key, value in data.items(): + full_key = f"{parent_key}.{key}" if parent_key else key + if isinstance(value, dict): # 处理字典 + form_html += f"

{key}

{render_form(value, full_key)}
" + else: + form_html += f"" + form_html += f"
" + return form_html + + +@app.route('/', methods=['GET', 'POST']) +def index(): + if request.method == 'POST': + data = load_yaml() + for key in request.form: + keys = key.split('.') + temp = data + for k in keys[:-1]: + temp = temp[k] + temp[keys[-1]] = request.form[key] + save_yaml(data) + return redirect('/') + + data = load_yaml() + form_html = render_form(data) + return render_template('index_overall.html', form_html=form_html) + + +if __name__ == '__main__': + app.run(debug=True) diff --git a/tools/setConfig/templates/config_overall.yaml b/tools/setConfig/templates/config_overall.yaml new file mode 100644 index 0000000000000000000000000000000000000000..d22dc6847dd3ee0ade13373233a91e38ec1f5770 --- /dev/null +++ b/tools/setConfig/templates/config_overall.yaml @@ -0,0 +1,145 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# Overall settings +# CSST-Sim Group, 2024/01/08 +# +############################################### + +# Base diretories and naming setup +# can add some of the command-line arguments here as well; +# ok to pass either way or both, as long as they are consistent +work_dir: "/public/home/fangyuedong/project/workplace/" +run_name: "ext_on" + +# Project cycle and run counter are used to name the outputs +project_cycle: 9 +run_counter: 1 + +# Run options +run_option: + # Output catalog only? + # If yes, no imaging simulation will be run. Only the catalogs + # of corresponding footprints will be generated. + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure the input catalog: options should be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/" + star_cat: "starcat_C9/" + galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" + + SED_templates_path: + star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/starcat_C9/" + galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/" + AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: NO + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + # Whether to apply milky way extinction to galaxies + enable_mw_ext_gal: YES + planck_ebv_map: "/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits" + +############################################### +# Observation setting +############################################### +obs_setting: + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat" + + obs_config_file: "/public/home/fangyuedong/project/csst_msc_sim/config/obs_config_SCI.yaml" + + # 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, 1, 2, 3, 4] + + # Whether to enable astrometric modeling + enable_astrometric_model: YES + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + # mag_sat_margin: -15. + + # 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 models for photometry survey simulation + psf_pho_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/dataC6/psfCube1" + # PSF models for slitless spectrum survey simulation + psf_sls_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SLS_PSF_PCA_fp/" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": get shear values from catalog + shear_type: "constant" + + # For constant shear field + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Output options +############################################### +output_setting: + output_format: "channels" # Whether to export as 16 channels (subimages) with pre- and over-scan ("image"/"channels") + shutter_output: NO # Whether to export shutter effect 16-bit image + prnu_output: NO # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + +############################################### +# 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/tools/setConfig/templates/index_obs.html b/tools/setConfig/templates/index_obs.html new file mode 100644 index 0000000000000000000000000000000000000000..5199558d997fceadbd7caa31daa7f39b19508254 --- /dev/null +++ b/tools/setConfig/templates/index_obs.html @@ -0,0 +1,246 @@ + + + + + + csst_msc_sim_CONF + + + + +
+

配置obs_config_SCI.yaml

+

说明: + 1,该脚本可用于生成CSST主巡天成像仿真的参数文件,相关参数说明见本页脚注;2,用户必须修改相关路径参数，其他参数可参考默认值 +. +

+
+
+ {{ form_html | safe }} + +
+
+ + + + + diff --git a/tools/setConfig/templates/index_overall.html b/tools/setConfig/templates/index_overall.html new file mode 100644 index 0000000000000000000000000000000000000000..d61ccf82cd4be5e981c9b57a09eb19bc9d64605d --- /dev/null +++ b/tools/setConfig/templates/index_overall.html @@ -0,0 +1,258 @@ + + + + + + csst_msc_sim_CONF + + + + +
+

配置config_overall.yaml

+

说明: + 1,该脚本可用于生成CSST主巡天成像仿真的参数文件,相关参数说明见本页脚注;2,用户必须修改相关路径参数，其他参数可参考默认值. +

+
+
+ {{ form_html | safe }} + +
+
+ + + + + diff --git a/tools/setConfig/templates/obs_config_SCI.yaml b/tools/setConfig/templates/obs_config_SCI.yaml new file mode 100644 index 0000000000000000000000000000000000000000..ef82bcdc6c0e17cc788ca396dd4fa3be087a159a --- /dev/null +++ b/tools/setConfig/templates/obs_config_SCI.yaml @@ -0,0 +1,85 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# For single exposure type: +# SCI-WIDE +# CSST-Sim Group, 2024/01/08 +# +############################################### + +# Observation type +obs_type: "SCI" +obs_type_code: "101" +obs_id: "00000001" # this setting will only be used if pointing list file is not given + +# Define list of chips +# run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips +#run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips +run_chips: [17, 22] + +# Define observation sequence +call_sequence: + # Accumulate fluxes from objects + scie_obs: + # [Optional]: exposure time of the pointing will be used as default. + # Set it here is you want to override the default + # exptime: 150. # [s] + shutter_effect: YES + flat_fielding: YES + field_dist: YES + # Accumulate fluxes from sky background + sky_background: + # [Optional]: exposure time of the pointing will be used as default. + # Set it here is you want to override the default + # exptime: 150. # [s] + shutter_effect: YES + flat_fielding: YES + enable_straylight_model: YES + # flat_level: set the total skybackground value (e-) in the exptime,if none,set null, or delete the key + # flat_level_filt: the vale of "flat_level" is in the filter "flat_level_filt", can set NUV, u, g, r, i, z, y, if + # none,set null,or delete the key + flat_level: null + flat_level_filt: null + # Apply PRNU to accumulated photons + PRNU_effect: {} + # Accumulate photons caused by cosmic rays + cosmic_rays: + # [Optional]: exposure time of the pointing will be used as default. + # Set it here is you want to override the default + # exptime: 150. # [s] + save_cosmic_img: YES # # Whether to export cosmic ray image + # Add Poission noise and dark current + poisson_and_dark: + # [Optional]: exposure time of the pointing will be used as default. + # Set it here is you want to override the default + # exptime: 150. # [s] + add_dark: YES + # Simulate brighter fatter effects + bright_fatter: {} + # Add detector defects: hot/warm pixels, bad columns + detector_defects: + hot_pixels: YES + dead_pixels: YES + bad_columns: YES + # Apply response nonlinearity + nonlinearity: {} + # Apply CCD Saturation & Blooming + blooming: {} + # Run CTE simulation + # CTE_effect: {} + # Add prescan and overscan + prescan_overscan: + add_dark: YES + # Add bias + bias: + bias_16channel: YES + # Add readout noise + readout_noise: {} + # Apply gain + gain: + gain_16channel: YES + # Output the final image + quantization_and_output: + format_output: YES +...