diff --git a/Catalog/C6_50sqdeg_ns.py b/Catalog/C6_50sqdeg_ns.py new file mode 100644 index 0000000000000000000000000000000000000000..3a80f2802d85cb5e74d29ae734eecad87b41dd04 --- /dev/null +++ b/Catalog/C6_50sqdeg_ns.py @@ -0,0 +1,526 @@ +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 ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar +from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist +from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position + +# (TEST) +from astropy.cosmology import FlatLambdaCDM +from astropy import constants +from astropy import units as U +from astropy.coordinates import SkyCoord +from astropy.io import fits +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 + +bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5'] + +qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits'] + +# star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5'] +star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)] + +star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5'] + +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_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): + assert NSIDE == 2**healpixOrder + shift = healpixOrder - bundleOrder + shift = 2*shift + + nside_bundle = 2**bundleOrder + nside_healpix= 2**healpixOrder + + healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) + bundleID_nest = (healpixID_nest >> shift) + bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) + + return bundleID_ring + +def get_agnsed_file(bundle_file_name): + return qsosed_file_list[bundle_file_list.index(bundle_file_name)] + +def get_star_cat(ra_pointing, dec_pointing): + pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg) + max_dist = 10 + return_star_path = None + for star_file, center in zip(star_file_list, star_center_list): + center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg) + dist = pointing_c.separation(center_c).to(U.deg).value + if dist < max_dist: + return_star_path = star_file + max_dist = dist + return return_star_path + +class Catalog(CatalogBase): + def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): + super().__init__() + self.cat_dir = 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. + + if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: + # Get the cloest star catalog file + star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec) + star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name) + 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, 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 + ) + + 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] + + 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.filter_type == 'NUV': + param['mag_use_normal'] = gals['mag_csst_nuv'][igals] + else: + param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] + if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): + continue + + param['z'] = gals['redshift'][igals] + param['model_tag'] = 'None' + param['g1'] = gals['shear'][igals][0] + param['g2'] = gals['shear'][igals][1] + param['kappa'] = gals['kappa'][igals] + param['e1'] = gals['ellipticity_true'][igals][0] + param['e2'] = gals['ellipticity_true'][igals][1] + + # For shape calculation + param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity( + e1=gals['ellipticity_true'][igals][0], + e2=gals['ellipticity_true'][igals][1], + rotation=self.rotation, + unit='radians') + # param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) + if param['ell_total'] > 0.9: + continue + # phi_e = cmath.phase(complex(param['e1'], param['e2'])) + # param['e1'] = param['ell_total'] * np.cos(phi_e + 2*self.rotation) + # param['e2'] = param['ell_total'] * np.sin(phi_e + 2*self.rotation) + + param['e1_disk'] = param['e1'] + param['e2_disk'] = param['e2'] + param['e1_bulge'] = param['e1'] + param['e2_bulge'] = param['e2'] + + + param['delta_ra'] = 0 + param['delta_dec'] = 0 + + # Masses + param['bulgemass'] = gals['bulgemass'][igals] + param['diskmass'] = gals['diskmass'][igals] + + param['size'] = gals['size'][igals] + if param['size'] > self.max_size: + self.max_size = param['size'] + + # Sersic index + param['disk_sersic_idx'] = 1. + param['bulge_sersic_idx'] = 4. + + # Sizes + param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) + if param['bfrac'] >= 0.6: + param['hlr_bulge'] = param['size'] + param['hlr_disk'] = param['size'] * (1. - param['bfrac']) + else: + param['hlr_disk'] = param['size'] + param['hlr_bulge'] = param['size'] * param['bfrac'] + + # SED coefficients + param['coeff'] = gals['coeff'][igals] + param['detA'] = gals['detA'][igals] + + # Others + param['galType'] = gals['type'][igals] + param['veldisp'] = gals['veldisp'][igals] + + # TEST no redening and no extinction + param['av'] = 0.0 + param['redden'] = 0 + + # 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_C6_bundle{:06}.h5".format(bundleID) + file_path = os.path.join(self.galaxy_path, bundle_file) + gals_cat = h5.File(file_path, 'r')['galaxies'] + gals = gals_cat[str(pix)] + + # Get corresponding AGN SED file + agnsed_file = get_agnsed_file(bundle_file) + agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file) + self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data + + self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file) + + del gals + except Exception as e: + traceback.print_exc() + self.logger.error(str(e)) + print(e) + + if self.logger is not None: + self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) + self.logger.info("number of objects in catalog: %d"%(len(self.objs))) + else: + print("number of objects in catalog: ", len(self.objs)) + + def load_sed(self, obj, **kwargs): + if obj.type == 'star': + # _, wave, flux = tag_sed( + # h5file=self.tempSED_star, + # model_tag=obj.param['model_tag'], + # teff=obj.param['teff'], + # logg=obj.param['logg'], + # feh=obj.param['feh'] + # ) + 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) + lamb = np.arange(2000, 11001+0.5, 0.5) + y = speci(lamb) + # erg/s/cm2/A --> photon/s/m2/A + all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 + sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) + + if obj.type == 'quasar': + # integrate to get the magnitudes + sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T + sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') + sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) + interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) + obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) + # mag = getABMAG(interFlux, self.filt.bandpass_full) + # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) + del wave + del flux + return sed diff --git a/Catalog/data/starSpecInterp.so b/Catalog/data/starSpecInterp.so new file mode 100755 index 0000000000000000000000000000000000000000..150d42c5d743164a22eef536f846317e85aae387 Binary files /dev/null and b/Catalog/data/starSpecInterp.so differ diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 8253567f3d9f1eee004846bca44d6e17a6e60006..f4f14351df9d94ff048dec8212ffcecf67c342af 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -19,7 +19,7 @@ from astropy.time import Time from astropy import wcs from ObservationSim.Config._util import get_obs_id, get_file_type -from datetime import datetime +from datetime import datetime, timezone # import socket import platform import toml @@ -352,7 +352,10 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointing_id = '00000001', po # ccdnum = str(k) datetime_obs = datetime.utcfromtimestamp(time_pt) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) + # print(datetime_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]) datetime_obs = datetime.utcfromtimestamp(np.round(datetime_obs.timestamp(), 1)) + # print(datetime_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]) # date_obs = datetime_obs.strftime("%y%m%d") # time_obs = datetime_obs.strftime("%H%M%S%f")[:-5] @@ -551,30 +554,32 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p h_ext['PIXSCAL1'] = pixel_scale h_ext['PIXSCAL2'] = pixel_scale h_ext['EXPTIME'] = exptime - h_ext['DARKTIME'] = exptime + readoutTime + h_ext['DARKTIME'] = exptime datetime_obs = datetime.utcfromtimestamp(timestamp) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) tstart = Time(datetime_obs) t_shutter_os = tstart t_shutter_oe = Time(tstart.mjd + t_shutter_open / 86400., format="mjd") t_shutter_co = Time(tstart.mjd + exptime / 86400., format="mjd") t_shutter_ce = Time(tstart.mjd + (exptime + t_shutter_close) / 86400., format="mjd") - t_shutter_os1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_os.unix).timestamp(), 1)) + t_shutter_os1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_os.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['SHTOPEN0'] = t_shutter_os1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] - t_shutter_oe1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_oe.unix).timestamp(), 1)) + t_shutter_oe1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_oe.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['SHTOPEN1'] = t_shutter_oe1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] - t_shutter_co1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_co.unix).timestamp(), 1)) + t_shutter_co1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_co.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['SHTCLOS0'] = t_shutter_co1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] - t_shutter_ce1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_ce.unix).timestamp(), 1)) + t_shutter_ce1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_ce.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['SHTCLOS1'] = t_shutter_ce1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] + tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd") tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd") # tstart1=tstart.datetime.replace(microsecond=round(tstart.datetime.microsecond, -5)) - tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart_read.unix).timestamp(), 1)) + tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart_read.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['ROTIME0'] = tstart1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] # tend_read1 = tend_read.datetime.replace(microsecond=round(tend_read.datetime.microsecond, -5)) - tend_read1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tend_read.unix).timestamp(), 1)) + tend_read1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tend_read.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) h_ext['ROTIME1'] = tend_read1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] # h_ext['POS_ANG'] = pa header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size, diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py index a5e7e6b84e2cb4a6970b6f0af01c4aad8884b386..a65943e7971b0c1ba8453b9783c4cc5df69578b8 100644 --- a/ObservationSim/MockObject/CatalogBase.py +++ b/ObservationSim/MockObject/CatalogBase.py @@ -56,6 +56,8 @@ class CatalogBase(metaclass=ABCMeta): "teff":0., "logg":0., "feh":0., + "DM":0., + "stellarMass":1., # C6 galaxies parameters "e1":0., "e2":0., diff --git a/ObservationSim/MockObject/FlatLED.py b/ObservationSim/MockObject/FlatLED.py index 6c9ee5e4f9ab1863ed5684bd5e756dae521a281f..21d0b688557bf5cf75ee32d540cc5a444c7902c5 100755 --- a/ObservationSim/MockObject/FlatLED.py +++ b/ObservationSim/MockObject/FlatLED.py @@ -11,6 +11,7 @@ from scipy.interpolate import griddata from astropy.table import Table from ObservationSim.MockObject.SpecDisperser import SpecDisperser from scipy import interpolate +import gc from ObservationSim.MockObject.MockObject import MockObject # from ObservationSim.Straylight import calculateSkyMap_split_g @@ -96,9 +97,48 @@ class FlatLED(object): U = griddata(X_, Z_, ( M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)], N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]), - method='cubic') + method='linear') U = U/np.mean(U) flatImage = U*fluxLED[led_type]*1000 + gc.collect() + return flatImage + + + ### + ### return LED flat, e/s + ### + def getLEDImage1(self, led_type='LED1'): + # cwave = cwaves[led_type] + flat = fits.open(os.path.join(self.flatDir, 'model_' + cwaves_name[led_type] + 'nm.fits')) + xlen = flat[0].header['NAXIS1'] + ylen = 601 + + i = self.chip.rowID - 1 + j = self.chip.colID - 1 + + x = np.linspace(0, self.chip.npix_x, int(xlen/6.)) + y = np.linspace(0, self.chip.npix_y, int(ylen/5.)) + xx, yy = np.meshgrid(x, y) + + a1 = flat[0].data[int(ylen*i/5.):int(ylen*i/5.)+int(ylen/5.), int(xlen*j/6.):int(xlen*j/6.)+int(xlen/6.)] + # z = np.sin((xx+yy+xx**2+yy**2)) + # fInterp = interp2d(xx, yy, z, kind='linear') + + X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None])) + Z_ = a1.flatten() + + n_x = np.arange(0, self.chip.npix_x , 1) + n_y = np.arange(0, self.chip.npix_y, 1) + + M, N = np.meshgrid(n_x, n_y) + + U = griddata(X_, Z_, ( + M[0:self.chip.npix_y, 0:self.chip.npix_x], + N[0:self.chip.npix_y, 0:self.chip.npix_x ]), + method='linear') + U = U/np.mean(U) + flatImage = U*fluxLED[led_type]*1000 + gc.collect() return flatImage def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]): @@ -114,7 +154,9 @@ class FlatLED(object): for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] - unitFlatImg = self.getLEDImage(led_type=led_type) + # unitFlatImg = self.getLEDImage(led_type=led_type) + unitFlatImg = self.getLEDImage1(led_type=led_type) + # print("---------------TEST mem:",np.mean(unitFlatImg)) led_wave = cwaves[led_type] led_fwhm = cwaves_fwhm[led_type] led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm) @@ -134,6 +176,7 @@ class FlatLED(object): ledStat = ledStat[0:int(led_type[3:])-1]+nledStat+ledStat[int(led_type[3:]):] ledTimes[int(led_type[3:])-1] = exp_t * 1000 + gc.collect() return ledFlat, ledStat, ledTimes def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]): @@ -150,7 +193,9 @@ class FlatLED(object): for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] - unitFlatImg = self.getLEDImage(led_type=led_type) + # unitFlatImg = self.getLEDImage(led_type=led_type) + unitFlatImg = self.getLEDImage1(led_type=led_type) + # print("---------------TEST mem:",np.mean(unitFlatImg)) ledFlat_ = unitFlatImg*exp_t ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type] ledFlat_.astype(np.float32) diff --git a/ObservationSim/sim_steps/add_LED_flat.py b/ObservationSim/sim_steps/add_LED_flat.py index bacfab4b9b4617bea33d5d335f46329dd4c75bda..751485d79cde657158e53c1dea3f580850b3b1d2 100644 --- a/ObservationSim/sim_steps/add_LED_flat.py +++ b/ObservationSim/sim_steps/add_LED_flat.py @@ -2,6 +2,11 @@ import numpy as np from ObservationSim.MockObject import FlatLED import galsim +from astropy.time import Time +from datetime import datetime, timezone + +import gc + def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param): if not hasattr(self, 'h_ext'): @@ -17,6 +22,7 @@ def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param): pf_map = led_flat self.updateHeaderInfo(header_flag='ext', keys = ['LEDSTAT'], values = [ledstat]) self.updateHeaderInfo(header_flag='ext', keys = ['LEDT01','LEDT02','LEDT03','LEDT04','LEDT05','LEDT06','LEDT07','LEDT08','LEDT09','LEDT10','LEDT11','LEDT12','LEDT13','LEDT14'], values = letts) + if obs_param["shutter_effect"] == True: pf_map = pf_map * chip.shutter_img @@ -26,4 +32,20 @@ def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param): self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) chip.img = chip.img + pf_map + + # renew header info + datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) + t_obs = Time(datetime_obs) + + ##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光 + t_obs_renew = Time(t_obs.mjd - (2.+0.5 +0.5) / 86400., format="mjd") + + t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) + self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) + + #dark time : 曝光时间+刷新后等带时间0.5s+点亮灯后0.5s+关闭快门时间1.5s+管快门后关灯前0.5+关灯后读出前等待0.5s + self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+0.5+1.5+0.5+0.5+pointing.exp_time]) + + gc.collect() return chip, filt, tel, pointing \ No newline at end of file diff --git a/ObservationSim/sim_steps/add_objects.py b/ObservationSim/sim_steps/add_objects.py index a6395c90adac68bb47deb8ed42208c206dd87b76..866d5389a82f9ebd7a4905e46316ce7f84dd13ec 100644 --- a/ObservationSim/sim_steps/add_objects.py +++ b/ObservationSim/sim_steps/add_objects.py @@ -7,6 +7,9 @@ import galsim from ObservationSim._util import get_shear_field from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS +from astropy.time import Time +from datetime import datetime, timezone + def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Get exposure time @@ -59,6 +62,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Get chip WCS if not hasattr(self, 'h_ext'): _, _ = self.prepare_headers(chip=chip, pointing=pointing) + chip_wcs = galsim.FitsWCS(header = self.h_ext) # Loop over objects @@ -211,4 +215,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): chip.img *= flat_normal del flat_normal + + # renew header info + datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) + t_obs = Time(datetime_obs) + + ##ccd刷新2s,等待0.5s，开始曝光 + t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd") + + t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) + self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) + + #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s + self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time]) + return chip, filt, tel, pointing \ No newline at end of file diff --git a/ObservationSim/sim_steps/add_sky_background.py b/ObservationSim/sim_steps/add_sky_background.py index 8b7d67e6f48d9e48aac46914676ef7a934231d8f..50a1731e8497e66e9503657a087d4b9f87c02610 100644 --- a/ObservationSim/sim_steps/add_sky_background.py +++ b/ObservationSim/sim_steps/add_sky_background.py @@ -3,6 +3,9 @@ import galsim from ObservationSim.Straylight import calculateSkyMap_split_g from ObservationSim.Instrument import FilterParam +from astropy.time import Time +from datetime import datetime, timezone + def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param): # Get exposure time @@ -120,4 +123,20 @@ def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param): chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param) else: chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param) + + # renew header info + datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) + t_obs = Time(datetime_obs) + + ##ccd刷新2s,等待0.5s，开始曝光 + t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd") + + t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) + self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) + + #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s + self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time]) + + return chip, filt, tel, pointing \ No newline at end of file diff --git a/ObservationSim/sim_steps/readout_output.py b/ObservationSim/sim_steps/readout_output.py index 89b034308f4f5caba74193d93bf37641cc7e19db..4147665de4c24f9ab540e512281c78ade87a0e70 100644 --- a/ObservationSim/sim_steps/readout_output.py +++ b/ObservationSim/sim_steps/readout_output.py @@ -5,6 +5,9 @@ from astropy.io import fits from ObservationSim.Instrument.Chip import ChipUtils as chip_utils from ObservationSim.Instrument.Chip import Effects +from astropy.time import Time +from datetime import datetime, timezone + def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param): self.chip_output.Log_info("Apply pre/over-scan") chip.img = chip_utils.AddPreScan(GSImage=chip.img, @@ -42,7 +45,17 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) if not hasattr(self, 'h_ext'): _, _ = self.prepare_headers(chip=chip, pointing=pointing) - self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1','EXPTIME'], values = [False,'','','','',0.0]) + # renew header info + datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) + t_obs = Time(datetime_obs) + + ##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光 + t_obs_renew = Time(t_obs.mjd - 2. / 86400., format="mjd") + + t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) + self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) gains1 = list(chip.gain_channel[0:8]) gains2 = list(chip.gain_channel[8:]) diff --git a/config/config_overall_newStar.yaml b/config/config_overall_newStar.yaml new file mode 100644 index 0000000000000000000000000000000000000000..10cdeee291310bef3acb4cb6f9bcc866808dd934 --- /dev/null +++ b/config/config_overall_newStar.yaml @@ -0,0 +1,140 @@ +--- +############################################### +# +# 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: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/" +run_name: "QSO_50sqdeg_test" + +# Project cycle and run counter are used to name the outputs +project_cycle: 9 +run_counter: 1 + +# Run options +run_option: + use_mpi: NO + + # 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: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/" + star_cat: "star_catalog/" + # galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" + + SED_templates_path: + star_SED: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/star_catalog/" + # galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/" + # AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/" + + # Only simulate stars? + star_only: YES + + # Only simulate galaxies? + galaxy_only: NO + +############################################### +# 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: "/nfsdata/share/CSSOSDataProductsSims/data/pointing_50_1_n.dat" + + obs_config_file: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/csst-simulation/config/obs_config_SCI_WIDE_phot.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: [1] + + # Whether to enable astrometric modeling + enable_astrometric_model: NO + + # 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: "/nfsdata/share/CSSOSDataProductsSims/data/psfCube1" + # PSF models for slitless spectrum survey simulation + psf_sls_dir: "/nfsdata/share/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": 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/setup.py b/setup.py index 1a902d5258d8989ccda0d0e23d4e4d2a7f6708b8..ade3b6961a51c7e913e2ed19c96c5d82521a8c47 100644 --- a/setup.py +++ b/setup.py @@ -80,7 +80,7 @@ setup(name='CSSTSim', 'ObservationSim.Instrument.data.throughputs': ['*.txt', '*.dat'], 'ObservationSim.Instrument.data.sls_conf': ['*.conf', '*.fits'], 'ObservationSim.Instrument.data.flatCube': ['*.fits'], - 'Catalog.data': ['*.fits'], + 'Catalog.data': ['*.fits','*.so'], 'ObservationSim.Config.Header':['*.header','*.lst'], 'ObservationSim.Straylight.data': ['*.dat'], 'ObservationSim.Straylight.data.sky': ['*.dat'],