diff --git a/Catalog/Calib_Catalog.py b/Catalog/Calib_Catalog.py new file mode 100644 index 0000000000000000000000000000000000000000..3987ef548097419e2ca2889739e8e2134211d9ee --- /dev/null +++ b/Catalog/Calib_Catalog.py @@ -0,0 +1,673 @@ +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 "CALIB_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"][ + "CALIB_cat"] and not config["catalog_options"]["star_only"]: + self.CALIB_cat_path = os.path.join(config["data_dir"], + config["catalog_options"]["input_path"]["CALIB_cat"]) + self.CALIB_SED_path = os.path.join(config["data_dir"], + config["catalog_options"]["SED_templates_path"]["CALIB_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])) + # 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] + 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['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] + # 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_calibObj(self): + + data = Table.read(self.CALIB_cat_path) + ra_arr = data['RA'] + dec_arr = data['DEC'] + + ngals = len(data) + + # 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'] = data['RA'][igals] + param['dec_orig'] = data['DEC'][igals] + param['mag_use_normal'] = data['MAG_g'][igals] + # if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): + # continue + + param['z'] = -99 + param['model_tag'] = 'None' + param['g1'] = 0 + param['g2'] = 0 + param['kappa'] = 0 + 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'] = 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] + # if param['size'] > self.max_size: + # self.max_size = param['size'] + + # Sersic index + param['disk_sersic_idx'] = data['SERSIC_N'][igals] + param['bulge_sersic_idx'] = 1. + param['hlr_bulge'] = data['RE'][igals] + param['hlr_disk'] = data['RE'][igals] + param['bfrac'] = 0 + # TEST no redening and no extinction + param['av'] = 0.0 + param['redden'] = 0 + + param['star'] = 4 + + # 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'] = data['SPEC_FN'][igals][0:-5] + + if param['star'] == 4: + 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(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 "CALIB_cat" in self.config["catalog_options"]["input_path"] and \ + self.config["catalog_options"]["input_path"][ + "CALIB_cat"] and not self.config["catalog_options"]["star_only"]: + try: + self._load_calibObj() + 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) + elif obj.type == 'calib': + data = Table.read(os.path.join(self.CALIB_SED_path,obj.id+'.fits')) + obj_w = data['WAVELENGTH'] + obj_f = data['FLUX'] + input_delt_w = np.min(obj_w[1:]-obj_w[0:-1]) + if input_delt_w > 0.5: + lamb = np.arange(2000, 11000 + 0.5, 0.5) + speci = interpolate.interp1d(obj_w, obj_f) + y1 = speci(lamb) + else: + lamb = obj_w + y1 = obj_f + # erg/s/cm2/A --> photon/s/m2/A + y1_phot = y1 * lamb / (cons.h.value * cons.c.value) * 1e-13 + sed = Table(np.array([lamb, y1_phot]).T, + names=('WAVELENGTH', 'FLUX')) + + 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) + return sed + 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/extension_header.header b/ObservationSim/Config/Header/extension_header.header index 35ef9f6851795a2ebc5141eafaab98e5a0e79794..19bf1fd6ce51187d0426a8d7c85cc9bb81b8c882 100644 --- a/ObservationSim/Config/Header/extension_header.header +++ b/ObservationSim/Config/Header/extension_header.header @@ -1 +1 @@ -XTENSION= 'IMAGE ' / extension type BITPIX = 16 / bits per data value NAXIS = 2 / number of data axes NAXIS1 = 9216 / length of first array axis NAXIS2 = 9232 / length of second array axis PCOUNT = 0 GCOUNT = 1 EXTNAME = 'SCI ' EXTVER = 1 BSCALE = 1 BZERO = 32768 BUNIT = 'ADU ' / physical unit of array values COMMENT ================================================================== COMMENT Detector information COMMENT ================================================================== CAMERA = 'MS' / camera of main survey DETSN = '12345678' / detector serial number DETNAME = 'CCD' / detector type DETTEMP1= 173.0 / detector temperature at EXPSTART(in Kelvin) DETTEMP2= 173.0 / detector temperature at EXPEND(in Kelvin) DETTEMP3= 173.0 / detector temperature at READT1(in Kelvin) DETSIZE = '9560x9264' / detector size DATASECT= '9216x9232' / data section PIXSCAL1= 0.074 / pixel scale for axis 1 PIXSCAL2= 0.074 / pixel scale for axis 2 PIXSIZE1= 10 / pixel size for axis 1 (in um) PIXSIZE2= 10 / pixel size for axis 2 (in um) COMMENT ================================================================== COMMENT CCD chip information COMMENT ================================================================== CHIPID = '08' / chip ID CHIPLAB = 'y-1' / chip label FILTER = 'y' / filter name NCHAN = 16 / number of readout channels PSCAN1 = 27 / horizontal prescan width, per readout channel PSCAN2 = 8 / vertical prescan width, per readout channel OSCAN1 = 16 / horizontal overscan width,per readout channel OSCAN2 = 16 / vertical overscan width,per readout channel COMMENT ================================================================== COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS COMMENT ================================================================== WCSAXES = 2 / number of World Coordinate System axes CRPIX1 = -10017.0 / x-coordinate of reference pixel CRPIX2 = 24876.0 / y-coordinate of reference pixel CRVAL1 = 62.228226 / first axis value at reference pixel CRVAL2 = -42.316932 / second axis value at reference pixel CTYPE1 = 'RA---TAN' / the coordinate type for the first axis CTYPE2 = 'DEC--TAN' / the coordinate type for the second axis CD1_1 = 1.88602083707394E-05 / partial of first axis coordinate w.r.t.x CD1_2 = 8.17455836176000E-06 / partial of first axis coordinate w.r.t.y CD2_1 = -8.1745583617600E-06 / partial of second axis coordinate w.r.t.x CD2_2 = 1.88602083707394E-05 / partial of second axis coordinate w.r.t.y OTHERS = '' / COMMENT ================================================================== COMMENT Readout information COMMENT ================================================================== GAINLVL = '01' / gain level GAIN01 = 1.1 / gain (channel 01) GAIN02 = 1.1 / gain (channel 02) GAIN03 = 1.1 / gain (channel 03) GAIN04 = 1.1 / gain (channel 04) GAIN05 = 1.1 / gain (channel 05) GAIN06 = 1.1 / gain (channel 06) GAIN07 = 1.1 / gain (channel 07) GAIN08 = 1.1 / gain (channel 08) GAIN09 = 1.1 / gain (channel 09) GAIN10 = 1.1 / gain (channel 10) GAIN11 = 1.1 / gain (channel 11) GAIN12 = 1.1 / gain (channel 12) GAIN13 = 1.1 / gain (channel 13) GAIN14 = 1.1 / gain (channel 14) GAIN15 = 1.1 / gain (channel 15) GAIN16 = 1.1 / gain (channel 16) RON01 = 5.0 / read noise (channel 01) RON02 = 5.0 / read noise (channel 02) RON03 = 5.0 / read noise (channel 03) RON04 = 5.0 / read noise (channel 04) RON05 = 5.0 / read noise (channel 05) RON06 = 5.0 / read noise (channel 06) RON07 = 5.0 / read noise (channel 07) RON08 = 5.0 / read noise (channel 08) RON09 = 5.0 / read noise (channel 09) RON10 = 5.0 / read noise (channel 10) RON11 = 5.0 / read noise (channel 11) RON12 = 5.0 / read noise (channel 12) RON13 = 5.0 / read noise (channel 13) RON14 = 5.0 / read noise (channel 14) RON15 = 5.0 / read noise (channel 15) RON16 = 5.0 / read noise (channel 16) READT0 = '2024-00-00T00:00:00'/ readout start time(UTC) READT1 = '2024-00-00T00:00:00'/ readout end time(UTC) ROSPEED = 10.0 / readout speed (in MHz) EXPTIME = 150.0 / exposure duration DARKTIME= 150.0 / dark current time COMMENT ================================================================== COMMENT Shutter information COMMENT ================================================================== SHTSTAT = T / shutter status SHTOPEN0= 0.0 / shutter open time (begin) SHTOPEN1= 0.0 / shutter open time (end) SHTCLOS0= 0.0 / shutter close time (begin) SHTCLOS1= 0.0 / shutter close time (end) COMMENT ================================================================== COMMENT LED information COMMENT ================================================================== LEDFLAG = 0 / main/backup LED LEDSTAT = '00000000000000' / LED status LEDEXPT = 0.0 / LED flash time (s) LEDTEMP = 173.0 / LED temperature (in K) COMMENT ================================================================== COMMENT Other information COMMENT ================================================================== CHECKSUM= '''abcde''' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = '''abcde''' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END +XTENSION= 'IMAGE ' / extension type BITPIX = 16 / bits per data value NAXIS = 2 / number of data axes NAXIS1 = 9216 / length of first array axis NAXIS2 = 9232 / length of second array axis PCOUNT = 0 GCOUNT = 1 EXTNAME = 'SCI ' EXTVER = 1 BSCALE = 1 BZERO = 32768 BUNIT = 'ADU ' / physical unit of array values COMMENT ========================================================================COMMENT Detector information COMMENT ========================================================================CAMERA = 'MS' / camera of main survey DETSN = '12345678' / detector serial number DETNAME = 'CCD' / detector type DETTEMP1= 173.0 / detector temperature at EXPSTART(in Kelvin) DETTEMP2= 173.0 / detector temperature at EXPEND(in Kelvin) DETTEMP3= 173.0 / detector temperature at READT1(in Kelvin) DETSIZE = '9560x9264' / detector size DATASECT= '9216x9232' / data section PIXSCAL1= 0.074 / pixel scale for axis 1 PIXSCAL2= 0.074 / pixel scale for axis 2 PIXSIZE1= 10 / pixel size for axis 1 (in um) PIXSIZE2= 10 / pixel size for axis 2 (in um) COMMENT ========================================================================COMMENT CCD chip information COMMENT ========================================================================CHIPID = '08' / chip ID CHIPLAB = 'y-1' / chip label FILTER = 'y' / filter name NCHAN = 16 / number of readout channels PSCAN1 = 27 / horizontal prescan width, per readout channel PSCAN2 = 8 / vertical prescan width, per readout channel OSCAN1 = 16 / horizontal overscan width,per readout channel OSCAN2 = 16 / vertical overscan width,per readout channel COMMENT ========================================================================COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS COMMENT ========================================================================WCSAXES = 2 / number of World Coordinate System axes CRPIX1 = -10017.0 / x-coordinate of reference pixel CRPIX2 = 24876.0 / y-coordinate of reference pixel CRVAL1 = 62.228226 / first axis value at reference pixel CRVAL2 = -42.316932 / second axis value at reference pixel CTYPE1 = 'RA---TAN' / the coordinate type for the first axis CTYPE2 = 'DEC--TAN' / the coordinate type for the second axis CD1_1 = 1.88602083707394E-05 / partial of first axis coordinate w.r.t.x CD1_2 = 8.17455836176000E-06 / partial of first axis coordinate w.r.t.y CD2_1 = -8.1745583617600E-06 / partial of second axis coordinate w.r.t.x CD2_2 = 1.88602083707394E-05 / partial of second axis coordinate w.r.t.y OTHERS = '' / COMMENT ========================================================================COMMENT Readout information COMMENT ========================================================================GAINLVL = '01' / gain level GAIN01 = 1.1 / gain (channel 01) GAIN02 = 1.1 / gain (channel 02) GAIN03 = 1.1 / gain (channel 03) GAIN04 = 1.1 / gain (channel 04) GAIN05 = 1.1 / gain (channel 05) GAIN06 = 1.1 / gain (channel 06) GAIN07 = 1.1 / gain (channel 07) GAIN08 = 1.1 / gain (channel 08) GAIN09 = 1.1 / gain (channel 09) GAIN10 = 1.1 / gain (channel 10) GAIN11 = 1.1 / gain (channel 11) GAIN12 = 1.1 / gain (channel 12) GAIN13 = 1.1 / gain (channel 13) GAIN14 = 1.1 / gain (channel 14) GAIN15 = 1.1 / gain (channel 15) GAIN16 = 1.1 / gain (channel 16) RON01 = 5.0 / read noise (channel 01) RON02 = 5.0 / read noise (channel 02) RON03 = 5.0 / read noise (channel 03) RON04 = 5.0 / read noise (channel 04) RON05 = 5.0 / read noise (channel 05) RON06 = 5.0 / read noise (channel 06) RON07 = 5.0 / read noise (channel 07) RON08 = 5.0 / read noise (channel 08) RON09 = 5.0 / read noise (channel 09) RON10 = 5.0 / read noise (channel 10) RON11 = 5.0 / read noise (channel 11) RON12 = 5.0 / read noise (channel 12) RON13 = 5.0 / read noise (channel 13) RON14 = 5.0 / read noise (channel 14) RON15 = 5.0 / read noise (channel 15) RON16 = 5.0 / read noise (channel 16) READT0 = '2024-00-00T00:00:00'/ readout start time(UTC) READT1 = '2024-00-00T00:00:00'/ readout end time(UTC) ROSPEED = 10.0 / readout speed (in MHz) EXPTIME = 150.0 / exposure duration DARKTIME= 150.0 / dark current time COMMENT ========================================================================COMMENT Shutter information COMMENT ========================================================================SHTSTAT = T / shutter status SHTOPEN0= 0.0 / shutter open time (begin) SHTOPEN1= 0.0 / shutter open time (end) SHTCLOS0= 0.0 / shutter close time (begin) SHTCLOS1= 0.0 / shutter close time (end) COMMENT ========================================================================COMMENT LED information COMMENT ========================================================================LEDFLAG = 0 / main/backup LED LEDSTAT = '00000000000000' / LED status LEDEXPT = 0.0 / LED flash time (s) LEDTEMP = 173.0 / LED temperature (in K) COMMENT ========================================================================COMMENT Other information COMMENT ========================================================================CHECKSUM= '''abcde''' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = '''abcde''' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END diff --git a/ObservationSim/Config/Header/global_header.header b/ObservationSim/Config/Header/global_header.header index 523f589a8629e874414f26d7379db93bb0fbc2ef..815c5437a6a2173412767cf8c4de42ed01417a2c 100644 --- a/ObservationSim/Config/Header/global_header.header +++ b/ObservationSim/Config/Header/global_header.header @@ -1 +1 @@ -SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NEXTEND = 1 / number of array dimensions GROUPS = F / ' ' DATE = '2021-03-04T09:30:00'/ the date on which this file was written FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name FILETYPE= 'SCIE ' / observation type TELESCOP= 'CSST ' / telescope used to acquire data INSTRUME= 'MSC ' / instrument used to acquire data RADECSYS= 'ICRS ' / reference coordinates system EQUINOX = 2000.0 / FITSCREA= 'C6' / FITS create software version COMMENT ================================================================== COMMENT Object information COMMENT ================================================================== OBJECT = '00000000' / object name TARGET = '+000000000000' / target name (hhmmss+ddmmss) OBSID = '00000000' / observation ID OBJ_RA = 62.228226 / R.A. of the object (degrees) OBJ_DEC = -42.316932 / declination of the object (degrees) COMMENT ================================================================== COMMENT Telescope information COMMENT ================================================================== REFFRAME= 'CSSTGSC-1.0' / guide star catalog version DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss) SATESWV = '0001' / software version in the satellite EXPSTART= 59130.5 / exposure start time (MJD) CABSTART= 59130.5 / (MJD) SUNANGL0= 50.0 / angle between sun and opt axis at CABSTART MOONANG0= 30.0 / angle between moon and opt axis at CABSTART TEL_ALT0= 20.0 / angle between opt axis and Elimb at CABSTART POS_ANG0= 20.0 / angle between y axis and NP at CABSTART POSI0_X = 0.0 / the orbital position in X at CABSTART POSI0_Y = 0.0 / the orbital position in Y at CABSTART POSI0_Z = 0.0 / the orbital position in Z at CABSTART VELO0_X = 0.0 / the orbital velocity in X at CABSTART VELO0_Y = 0.0 / the orbital velocity in Y at CABSTART VELO0_Z = 0.0 / the orbital velocity in Z at CABSTART EULER0_1= 0.0 / euler angle 1 at CABSTART EULER0_2= 0.0 / euler angle 2 at CABSTART EULER0_3= 0.0 / euler angle 3 at CABSTART RA_PNT0 = 0.0 / RA of the pointing (degrees) at CABSTART DEC_PNT0= 0.0 / DEC of the pointing (degrees) at CABSTART EXPEND = 0.0 / exposure end time (MJD) CABEND = 0.0 / (MJD) SUNANGL1= 50.0 / angle between sun and opt axis at CABEND MOONANG1= 30.0 / angle between moon and opt axis at CABEND TEL_ALT1= 20.0 / angle between opt axis and Elimb at CABEND POS_ANG1= 20.0 / angle between y axis and NP at CABEND POSI1_X = 0.0 / the orbital position in X at CABEND POSI1_Y = 0.0 / the orbital position in Y at CABEND POSI1_Z = 0.0 / the orbital position in Z at CABEND VELO1_X = 0.0 / the orbital velocity in X at CABEND VELO1_Y = 0.0 / the orbital velocity in Y at CABEND VELO1_Z = 0.0 / the orbital velocity in Z at CABEND EULER1_1= 0.0 / euler angle 1 at CABEND EULER1_2= 0.0 / euler angle 2 at CABEND EULER1_3= 0.0 / euler angle 3 at CABEND RA_PNT1 = 0.0 / RA of the pointing (degrees) at CABEND DEC_PNT1= 0.0 / DEC of the pointing (degrees) at CABEND EXPTIME = 150.0 / exposure duration EPOCH = 2000.0 / coordinate epoch COMMENT Other information COMMENT ================================================================== CHECKSUM= 'abcdefg ' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = 'abcdefg ' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END +SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NEXTEND = 1 / number of array dimensions GROUPS = F / ' ' DATE = '2021-03-04T09:30:00'/ the date on which this file was written FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name FILETYPE= 'SCIE ' / observation type TELESCOP= 'CSST ' / telescope used to acquire data INSTRUME= 'MSC ' / instrument used to acquire data RADECSYS= 'ICRS ' / reference coordinates system EQUINOX = 2000.0 / FITSCREA= 'C6' / FITS create software version COMMENT ========================================================================COMMENT Object information COMMENT ========================================================================OBJECT = '00000000' / object name TARGET = '+000000000000' / target name (hhmmss+ddmmss) OBSID = '00000000' / observation ID OBJ_RA = 62.228226 / R.A. of the object (degrees) OBJ_DEC = -42.316932 / declination of the object (degrees) COMMENT ========================================================================COMMENT Telescope information COMMENT ========================================================================REFFRAME= 'CSSTGSC-1.0' / guide star catalog version DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss) SATESWV = '0001' / software version in the satellite EXPSTART= 59130.5 / exposure start time (MJD) CABSTART= 59130.5 / (MJD) SUNANGL0= 50.0 / angle between sun and opt axis at CABSTART MOONANG0= 30.0 / angle between moon and opt axis at CABSTART TEL_ALT0= 20.0 / angle between opt axis and Elimb at CABSTART POS_ANG0= 20.0 / angle between y axis and NP at CABSTART POSI0_X = 0.0 / the orbital position in X at CABSTART POSI0_Y = 0.0 / the orbital position in Y at CABSTART POSI0_Z = 0.0 / the orbital position in Z at CABSTART VELO0_X = 0.0 / the orbital velocity in X at CABSTART VELO0_Y = 0.0 / the orbital velocity in Y at CABSTART VELO0_Z = 0.0 / the orbital velocity in Z at CABSTART EULER0_1= 0.0 / euler angle 1 at CABSTART EULER0_2= 0.0 / euler angle 2 at CABSTART EULER0_3= 0.0 / euler angle 3 at CABSTART RA_PNT0 = 0.0 / RA of the pointing (degrees) at CABSTART DEC_PNT0= 0.0 / DEC of the pointing (degrees) at CABSTART EXPEND = 0.0 / exposure end time (MJD) CABEND = 0.0 / (MJD) SUNANGL1= 50.0 / angle between sun and opt axis at CABEND MOONANG1= 30.0 / angle between moon and opt axis at CABEND TEL_ALT1= 20.0 / angle between opt axis and Elimb at CABEND POS_ANG1= 20.0 / angle between y axis and NP at CABEND POSI1_X = 0.0 / the orbital position in X at CABEND POSI1_Y = 0.0 / the orbital position in Y at CABEND POSI1_Z = 0.0 / the orbital position in Z at CABEND VELO1_X = 0.0 / the orbital velocity in X at CABEND VELO1_Y = 0.0 / the orbital velocity in Y at CABEND VELO1_Z = 0.0 / the orbital velocity in Z at CABEND EULER1_1= 0.0 / euler angle 1 at CABEND EULER1_2= 0.0 / euler angle 2 at CABEND EULER1_3= 0.0 / euler angle 3 at CABEND RA_PNT1 = 0.0 / RA of the pointing (degrees) at CABEND DEC_PNT1= 0.0 / DEC of the pointing (degrees) at CABEND EXPTIME = 150.0 / exposure duration EPOCH = 2000.0 / coordinate epoch COMMENT Other information COMMENT ========================================================================CHECKSUM= 'abcdefg ' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = 'abcdefg ' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 81c29cd63ff9e0109194e072bba9cd24d6fe00c2..2bb552a42aaee5b8a7d89941481aace30e5bd1b9 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -275,8 +275,13 @@ class Galaxy(MockObject): bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp) bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge = bulge.shear(bulge_shape) - - gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk + + if self.bfrac == 0: + gal = disk + elif self.bfrac == 1: + gal = bulge + else: + gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk # (TEST) Random knots # knots = galsim.RandomKnots(npoints=100, profile=disk) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 8e361b9b2a79250f2bc05fc0c162382be581be6f..552a186a302e2cba8798668d1a739f9caa4be307 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -26,6 +26,11 @@ class MockObject(object): elif self.param["star"] == 3: self.type = "stamp" ###mock_stamp_END + ###for calibration + elif self.param["star"] == 4: + self.type = "calib" + ###END + self.sed = None self.fd_shear = None # Place holder for outputs diff --git a/config/config_C6_dev_calib.yaml b/config/config_C6_dev_calib.yaml new file mode 100644 index 0000000000000000000000000000000000000000..4075d72123f1bbbd17a293ca4993c0a7281b537e --- /dev/null +++ b/config/config_C6_dev_calib.yaml @@ -0,0 +1,236 @@ +--- +############################################### +# +# 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: "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/" +data_dir: "/Volumes/EAGET/C6_data/inputData/" +run_name: "C6_new_sim_2sq_run2" +project_cycle: 6 +run_counter: 1 + +# 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: + galaxy_cat: + AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + CALIB_cat: "calibration_data/GP/calibrationCat_CHIP1_GI.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" + CALIB_SED: "calibration_data/GP/GPTable/" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: NO + + # 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 + # "CALIBRATION": falt, bias, dark with or without postflash + survey_type: "Spectroscopic" + #"LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null + #'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', + #'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365' + LED_TYPE: ['LED5','LED3'] + LED_TIME: [1.,0.1] + # unit e- ,flat level + FLAT_LEVEL: 20000 + FLAT_LEVEL_FIL: 'g' + + # 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: "/Volumes/EAGET/C6_data/inputData/" + 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: [1] + + # 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 + 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: "Gauss" + + # 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" + psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" +############################################### +# 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: NO # 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 + add_prescan: OFF + 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: NO # Whether to add cosmic-ray + cray_differ: NO # 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 + format_output: OFF + + # 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 +...