diff --git a/Catalog/wcs_test_C6.py b/Catalog/wcs_test_C6.py new file mode 100644 index 0000000000000000000000000000000000000000..9300e3806ab80bc1be3f256d5b2c55ac518734d2 --- /dev/null +++ b/Catalog/wcs_test_C6.py @@ -0,0 +1,515 @@ +import os +import galsim +import random +import numpy as np +import h5py as h5 +import healpy as hp +import astropy.constants as cons +import traceback +from astropy.coordinates import spherical_to_cartesian +from astropy.table import Table +from scipy import interpolate +from datetime import datetime + +from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar +from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist +from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position + +# (TEST) +from astropy.cosmology import FlatLambdaCDM +from astropy import constants +from astropy import units as U + +try: + import importlib.resources as pkg_resources +except ImportError: + # Try backported to PY<37 'importlib_resources' + import importlib_resources as pkg_resources + +NSIDE = 128 + +def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): + assert NSIDE == 2**healpixOrder + shift = healpixOrder - bundleOrder + shift = 2*shift + + nside_bundle = 2**bundleOrder + nside_healpix= 2**healpixOrder + + healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) + bundleID_nest = (healpixID_nest >> shift) + bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) + + return bundleID_ring + +class Catalog(CatalogBase): + def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): + super().__init__() + self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) + self.seed_Av = config["catalog_options"]["seed_Av"] + + self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) + + self.chip_output = chip_output + self.filt = filt + self.logger = chip_output.logger + + with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: + self.normF_star = Table.read(str(filter_path)) + + self.config = config + self.chip = chip + self.pointing = pointing + + self.max_size = 0. + + if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: + star_file = config["catalog_options"]["input_path"]["star_cat"] + star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"] + self.star_path = os.path.join(self.cat_dir, star_file) + self.star_SED_path = os.path.join(config["data_dir"], star_SED_file) + self._load_SED_lib_star() + if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: + galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] + self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) + self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) + self._load_SED_lib_gals() + + if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: + AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] + self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) + self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) + self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) + self._load_SED_lib_AGN() + + if "rotateEll" in config["catalog_options"]: + self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) + 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, self.rotation, logger=self.logger) + + # Need to deal with additional output columns + obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., + param['bulgemass'], param['diskmass'], param['detA'], + param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], + param['galType'], param['veldisp']) + + self.objs.append(obj) + + def _load_stars(self, stars, pix_id=None): + nstars = len(stars['sourceID']) + # Apply astrometric modeling + ra_arr = stars["RA"][:] + dec_arr = stars["Dec"][:] + pmra_arr = stars['pmra'][:] + pmdec_arr = stars['pmdec'][:] + rv_arr = stars['RV'][:] + parallax_arr = stars['parallax'][:] + if self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = pmra_arr.tolist() + pmdec_list = pmdec_arr.tolist() + rv_list = rv_arr.tolist() + parallax_list = parallax_arr.tolist() + dt = datetime.utcfromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=nstars, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2000", + input_date_str=date_str, + input_time_str=time_str + ) + for istars in range(nstars): + # # (TEST) + # if istars > 100: + # break + + param = self.initialize_param() + param['ra'] = ra_arr[istars] + param['dec'] = dec_arr[istars] + param['ra_orig'] = stars["RA"][istars] + param['dec_orig'] = stars["Dec"][istars] + param['pmra'] = pmra_arr[istars] + param['pmdec'] = pmdec_arr[istars] + param['rv'] = rv_arr[istars] + param['parallax'] = parallax_arr[istars] + if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): + continue + # param['mag_use_normal'] = stars['app_sdss_g'][istars] + param['mag_use_normal'] = 20. + # if param['mag_use_normal'] >= 26.5: + # continue + self.ids += 1 + param['id'] = stars['sourceID'][istars] + param['sed_type'] = stars['sourceID'][istars] + param['model_tag'] = stars['model_tag'][istars] + param['teff'] = stars['teff'][istars] + param['logg'] = stars['grav'][istars] + param['feh'] = stars['feh'][istars] + param['z'] = 0.0 + param['star'] = 1 # Star + obj = Star(param, logger=self.logger) + + # Append additional output columns to the .cat file + obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], + 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) + + self.objs.append(obj) + + def _load_AGNs(self): + data = Table.read(self.AGN_path) + ra_arr = data['ra'] + dec_arr = data['dec'] + nAGNs = len(data) + if self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = np.zeros(nAGNs).tolist() + pmdec_list = np.zeros(nAGNs).tolist() + rv_list = np.zeros(nAGNs).tolist() + parallax_list = [1e-9] * nAGNs + dt = datetime.utcfromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=nAGNs, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2000", + input_date_str=date_str, + input_time_str=time_str + ) + for iAGNs in range(nAGNs): + param = self.initialize_param() + param['ra'] = ra_arr[iAGNs] + param['dec'] = dec_arr[iAGNs] + param['ra_orig'] = data['ra'][iAGNs] + param['dec_orig'] = data['dec'][iAGNs] + param['z'] = data['z'][iAGNs] + param['appMag'] = data['appMag'][iAGNs] + param['absMag'] = data['absMag'][iAGNs] + + # NOTE: this cut cannot be put before the SED type has been assigned + if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): + continue + + # TEST no redening and no extinction + param['av'] = 0.0 + param['redden'] = 0 + + param['star'] = 2 # Quasar + param['id'] = data['igmlos'][iAGNs] + + if param['star'] == 2: + obj = Quasar(param, logger=self.logger) + + # Append additional output columns to the .cat file + obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) + self.objs.append(obj) + + def _load(self, **kwargs): + self.objs = [] + self.ids = 0 + + if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: + star_cat = h5.File(self.star_path, 'r')['catalog'] + for pix in self.pix_list: + try: + stars = star_cat[str(pix)] + self._load_stars(stars, pix_id=pix) + del stars + except Exception as e: + self.logger.error(str(e)) + print(e) + + if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: + for pix in self.pix_list: + try: + bundleID = get_bundleIndex(pix) + file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) + gals_cat = h5.File(file_path, 'r')['galaxies'] + gals = gals_cat[str(pix)] + self._load_gals(gals, pix_id=pix, cat_id=bundleID) + del gals + except Exception as e: + traceback.print_exc() + self.logger.error(str(e)) + print(e) + + if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: + try: + self._load_AGNs() + except Exception as e: + traceback.print_exc() + self.logger.error(str(e)) + print(e) + + if self.logger is not None: + self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) + self.logger.info("number of objects in catalog: %d"%(len(self.objs))) + else: + print("number of objects in catalog: ", len(self.objs)) + + def load_sed(self, obj, **kwargs): + if obj.type == 'star': + _, wave, flux = tag_sed( + h5file=self.tempSED_star, + model_tag=obj.param['model_tag'], + teff=obj.param['teff'], + logg=obj.param['logg'], + feh=obj.param['feh'] + ) + elif obj.type == 'galaxy' or obj.type == 'quasar': + factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) + if obj.type == 'galaxy': + flux = np.matmul(self.pcs, obj.coeff) * factor + # if np.any(flux < 0): + # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) + flux[flux < 0] = 0. + sedcat = np.vstack((self.lamb_gal, flux)).T + sed_data = getObservedSED( + sedCat=sedcat, + redshift=obj.z, + av=obj.param["av"], + redden=obj.param["redden"] + ) + wave, flux = sed_data[0], sed_data[1] + elif obj.type == 'quasar': + flux = self.SED_AGN[int(obj.id)] * 1e-17 + # if np.any(flux < 0): + # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) + flux[flux < 0] = 0. + # sedcat = np.vstack((self.lamb_AGN, flux)).T + wave = self.lamb_AGN + # print("sed (erg/s/cm2/A) = ", sed_data) + # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) + else: + raise ValueError("Object type not known") + speci = interpolate.interp1d(wave, flux) + lamb = np.arange(2000, 11001+0.5, 0.5) + y = speci(lamb) + # erg/s/cm2/A --> photon/s/m2/A + all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 + sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) + + if obj.type == 'quasar': + # integrate to get the magnitudes + sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T + sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') + sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) + interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) + obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) + # if obj.param['mag_use_normal'] >= 30: + # print("obj ID = %d"%obj.id) + # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) + # print("integrated flux = %.7f"%(interFlux)) + # print("app mag = %.3f"%obj.param['appMag']) + # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) + # print("obj ID = %d"%obj.id) + # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) + # print("integrated flux = %.7f"%(interFlux)) + # print("app mag = %.3f"%obj.param['appMag']) + # print("abs mag = %.3f"%obj.param['absMag']) + # mag = getABMAG(interFlux, self.filt.bandpass_full) + # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) + del wave + del flux + return sed diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 8a9eafe0210f74bf5825a32f71a39b23f8c40879..8def8abef9762caf3ac9b495e7f09b95576b3fa6 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -287,7 +287,8 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r r_dat['CRPIX1'] = -x1 r_dat['CRPIX2'] = -y1 - cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[1] + # cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[1] + cd = np.array([[ pixel_scale, 0], [0, -pixel_scale]])/3600. cd_rot = rotate_CD_matrix(cd, pa_aper) r_dat['CD1_1'] = cd_rot[0,0] r_dat['CD1_2'] = cd_rot[0,1] diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 32b95579c1657bf3783685a01c93dd524f845353..f490506335a04fc07fa33c3f8e2c9b0eea40a4f0 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -115,7 +115,7 @@ class Chip(FocalPlane): self._getCRdata() # Define the sensor model - if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric": + if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric": self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func) else: self.sensor = galsim.Sensor() diff --git a/ObservationSim/Instrument/FocalPlane.py b/ObservationSim/Instrument/FocalPlane.py index 55ef478911366aefdc433a6cae340206eb9ff614..ddca78f7aae07b6c24d914d3b36ba7dc58c37bd8 100755 --- a/ObservationSim/Instrument/FocalPlane.py +++ b/ObservationSim/Instrument/FocalPlane.py @@ -86,15 +86,16 @@ class FocalPlane(object): if (xcen == None) or (ycen == None): xcen = self.cen_pix_x ycen = self.cen_pix_y + # dudx = -np.cos(img_rot.rad) * pix_scale + # dudy = -np.sin(img_rot.rad) * pix_scale + # dvdx = -np.sin(img_rot.rad) * pix_scale + # dvdy = +np.cos(img_rot.rad) * pix_scale + dudx = -np.cos(img_rot.rad) * pix_scale - dudy = -np.sin(img_rot.rad) * pix_scale + dudy = +np.sin(img_rot.rad) * pix_scale dvdx = -np.sin(img_rot.rad) * pix_scale - dvdy = +np.cos(img_rot.rad) * pix_scale + dvdy = -np.cos(img_rot.rad) * pix_scale - # dudx = +np.sin(img_rot.rad) * pix_scale - # dudy = +np.cos(img_rot.rad) * pix_scale - # dvdx = -np.cos(img_rot.rad) * pix_scale - # dvdy = +np.sin(img_rot.rad) * pix_scale moscen = galsim.PositionD(x=xcen, y=ycen) sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees) affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen) diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 3db13dcdbec118746a97e918e488bd0b0dbbd3be..7874d10b13b6b0414323b2d7fd72aa66fac2bd44 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -51,7 +51,8 @@ class Galaxy(MockObject): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return -1 for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -59,7 +60,8 @@ class Galaxy(MockObject): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return -1 ratio = sub/full @@ -78,12 +80,15 @@ class Galaxy(MockObject): gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk gal = gal.withFlux(nphotons) + if fd_shear is not None: + g1 += fd_shear.g1 + g2 += fd_shear.g2 gal_shear = galsim.Shear(g1=g1, g2=g2) gal = gal.shear(gal_shear) gal = galsim.Convolve(psf, gal) - if fd_shear is not None: - gal = gal.shear(fd_shear) + # if fd_shear is not None: + # gal = gal.shear(fd_shear) objs.append(gal) final = galsim.Sum(objs) return final @@ -97,7 +102,8 @@ class Galaxy(MockObject): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return 2, None nphotons_sum = 0 @@ -131,6 +137,13 @@ class Galaxy(MockObject): real_wcs_local = self.real_wcs.local(self.real_pos) + disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp) + disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) + disk = disk.shear(disk_shape) + bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp) + bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) + bulge = bulge.shear(bulge_shape) + for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -138,7 +151,8 @@ class Galaxy(MockObject): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) # return False continue @@ -153,40 +167,49 @@ class Galaxy(MockObject): # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) - disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp) - disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) - disk = disk.shear(disk_shape) - bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp) - bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) - bulge = bulge.shear(bulge_shape) + # disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp) + # disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) + # disk = disk.shear(disk_shape) + # bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp) + # bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) + # bulge = bulge.shear(bulge_shape) + + gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk - gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk + gal_temp = gal_temp.withFlux(nphotons) + if not big_galaxy: # Not apply PSF for very big galaxy + gal_temp = galsim.Convolve(psf, gal_temp) + + if i == 0: + gal = gal_temp + else: + gal = gal + gal_temp # (TEST) Random knots # knots = galsim.RandomKnots(npoints=100, profile=disk) # kfrac = np.random.random()*(1.0 - self.bfrac) # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots - gal = gal.withFlux(nphotons) - gal_shear = galsim.Shear(g1=g1, g2=g2) - gal = gal.shear(gal_shear) + # gal = gal.withFlux(nphotons) + # gal_shear = galsim.Shear(g1=g1, g2=g2) + # gal = gal.shear(gal_shear) - if not big_galaxy: # Not apply PSF for very big galaxy - gal = galsim.Convolve(psf, gal) - if fd_shear is not None: - gal = gal.shear(fd_shear) + # if not big_galaxy: # Not apply PSF for very big galaxy + # gal = galsim.Convolve(psf, gal) + # # if fd_shear is not None: + # # gal = gal.shear(fd_shear) - # Use (explicit) stamps to draw - stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) + # # Use (explicit) stamps to draw + # stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) - xmax = max(xmax, stamp.xmax - stamp.xmin) - ymax = max(ymax, stamp.ymax - stamp.ymin) + # xmax = max(xmax, stamp.xmax - stamp.xmin) + # ymax = max(ymax, stamp.ymax - stamp.ymin) - photons = stamp.photons - photons.x += x_nominal - photons.y += y_nominal - photons_list.append(photons) - del gal + # photons = stamp.photons + # photons.x += x_nominal + # photons.y += y_nominal + # photons_list.append(photons) + # del gal # # [C6 TEST] # print('xmax = %d, ymax = %d '%(xmax, ymax)) @@ -195,8 +218,21 @@ class Galaxy(MockObject): # top_stats = snapshot.statistics('lineno') # for stat in top_stats[:10]: # print(stat) + if fd_shear: + g1 += fd_shear.g1 + g2 += fd_shear.g2 + gal_shear = galsim.Shear(g1=g1, g2=g2) + gal = gal.shear(gal_shear) + # if fd_shear is not None: + # gal = gal.shear(fd_shear) + + stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) + photons = stamp.photons + photons.x += x_nominal + photons.y += y_nominal + photons_list.append(photons) - stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) + # stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) stamp.wcs = real_wcs_local stamp.setCenter(x_nominal, y_nominal) bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) @@ -226,7 +262,8 @@ class Galaxy(MockObject): else: # Return code 0: object photons missed this detector print("obj %s missed"%(self.id)) - self.logger.info("obj %s missed"%(self.id)) + if self.logger: + self.logger.info("obj %s missed"%(self.id)) return 0, pos_shear # # [C6 TEST] @@ -297,14 +334,17 @@ class Galaxy(MockObject): # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots gal = gal.withFlux(tel.pupil_area * exptime) + if fd_shear: + g1 += fd_shear.g1 + g2 += fd_shear.g2 gal_shear = galsim.Shear(g1=g1, g2=g2) gal = gal.shear(gal_shear) gal = galsim.Convolve(psf, gal) if not big_galaxy: # Not apply PSF for very big galaxy gal = galsim.Convolve(psf, gal) - if fd_shear is not None: - gal = gal.shear(fd_shear) + # if fd_shear is not None: + # gal = gal.shear(fd_shear) starImg = gal.drawImage(wcs=real_wcs_local, offset=offset) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 0b5ff7fc269577841ac64a043b7481ef49001644..102191e15af187ca7566edd6b45654ca3df875e3 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -32,6 +32,7 @@ class MockObject(object): # Place holder for outputs self.additional_output_str = "" + self.fd_shear = None self.logger = logger @@ -61,7 +62,7 @@ class MockObject(object): dec = self.param["dec"] return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees) - def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, img_header=None): + def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None): self.posImg = img.wcs.toImage(self.getPosWorld()) self.localWCS = img.wcs.local(self.posImg) if (fdmodel is not None) and (chip is not None): @@ -79,8 +80,10 @@ class MockObject(object): dx = x - self.x_nominal dy = y - self.y_nominal self.offset = galsim.PositionD(dx, dy) - - if img_header is not None: + + if chip_wcs is not None: + self.real_wcs = chip_wcs + elif img_header is not None: self.real_wcs = galsim.FitsWCS(header=img_header) else: self.real_wcs = None @@ -132,7 +135,8 @@ class MockObject(object): full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) return 2, None nphotons_sum = 0 @@ -149,6 +153,8 @@ class MockObject(object): self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, img_real_wcs=self.real_wcs) + print(self.real_pos.x, self.real_pos.y) + x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) y_nominal = int(np.floor(y + 0.5)) @@ -164,7 +170,8 @@ class MockObject(object): sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) + if self.logger: + self.logger.error(e) # return False continue @@ -214,7 +221,8 @@ class MockObject(object): else: # Return code 0: object photons missed this detector print("obj %s missed"%(self.id)) - self.logger.info("obj %s missed"%(self.id)) + if self.logger: + self.logger.info("obj %s missed"%(self.id)) return 0, pos_shear del photons_list diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 3c9af7863dd84ac2370b49c8d7a20802409f9799..785370b5cdf11e0aca2870675065612a940534ef 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -152,7 +152,7 @@ class Observation(object): cut_filter = temp_filter if self.config["ins_effects"]["field_dist"] == True: - self.fd_model = FieldDistortion(chip=chip) + self.fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) else: self.fd_model = None @@ -177,6 +177,8 @@ class Observation(object): xcen=chip.x_cen, ycen=chip.y_cen, extName='raw') + + chip_wcs = galsim.FitsWCS(header=h_ext) for j in range(self.nobj): @@ -208,7 +210,7 @@ class Observation(object): continue # [TODO] Testing - # chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()])) + chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()])) # Exclude very bright/dim objects (for now) if cut_filter.is_too_bright( @@ -245,7 +247,9 @@ class Observation(object): raise ValueError("Unknown shear input") # Get position of object on the focal plane - pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=h_ext) + pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=h_ext) + + print(pos_img.x, pos_img.y) # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # Otherwise they will be considered missed objects diff --git a/ObservationSim/PSF/FieldDistortion.py b/ObservationSim/PSF/FieldDistortion.py index 95dc40e7547e57e8e1545c42ac6444708271e09b..f7ee0f8fc381c8a1cb14488aab56bc6dd4acf6b7 100644 --- a/ObservationSim/PSF/FieldDistortion.py +++ b/ObservationSim/PSF/FieldDistortion.py @@ -1,9 +1,10 @@ import galsim import numpy as np +import cmath class FieldDistortion(object): - def __init__(self, chip, fdModel=None, fdModel_path=None): + def __init__(self, chip, fdModel=None, fdModel_path=None, img_rot=0.): if fdModel is None: if hasattr(chip, 'fdModel'): self.fdModel = chip.fdModel @@ -15,6 +16,7 @@ class FieldDistortion(object): raise ValueError("Error: no field distortion model has been specified!") else: self.fdModel = fdModel + self.img_rot = img_rot self.ifdModel = self.fdModel["wave1"] self.ixfdModel = self.ifdModel["xImagePos"] self.iyfdModel = self.ifdModel["yImagePos"] @@ -42,7 +44,7 @@ class FieldDistortion(object): return False return True - def get_distorted(self, chip, pos_img, bandpass=None): + def get_distorted(self, chip, pos_img, bandpass=None, img_rot=None): """ Get the distored position for an undistorted image position Parameters: @@ -58,14 +60,14 @@ class FieldDistortion(object): """ if not self.isContainObj_FD(chip=chip, pos_img=pos_img): return galsim.PositionD(-1, -1), None + if not img_rot: + img_rot = np.radians(self.img_rot) + else: + img_rot = np.radians(img_rot) x, y = pos_img.x, pos_img.y x = self.ixfdModel(x, y)[0][0] y = self.iyfdModel(x, y)[0][0] - ix_dx = self.ifx_dx(x, y) - ix_dy = self.ifx_dy(x, y) - iy_dx = self.ify_dx(x, y) - iy_dy = self.ify_dy(x, y) - if self.irsModel is not None: + if self.irsModel: # x1LowI, x1UpI, y1LowI, y1UpI = self.irsModel["interpLimit"] # if (x1LowI-x)*(x1UpI-x) <=0 and (y1LowI-y)*(y1UpI-y)<=0: # dx = self.ixrsModel(x, y)[0][0] @@ -88,8 +90,25 @@ class FieldDistortion(object): ix_dy = self.ifx_dy(x, y) + self.irx_dy(x, y) iy_dx = self.ify_dx(x, y) + self.iry_dx(x, y) iy_dy = self.ify_dy(x, y) + self.iry_dy(x, y) + else: + ix_dx = self.ifx_dx(x, y) + ix_dy = self.ifx_dy(x, y) + iy_dx = self.ify_dx(x, y) + iy_dy = self.ify_dy(x, y) g1k_fd = 0.0 + (iy_dy - ix_dx) / (iy_dy + ix_dx) g2k_fd = 0.0 - (iy_dx + ix_dy) / (iy_dy + ix_dx) + + # [TODO] [TESTING] Rotate the shear: + g_abs = np.sqrt(g1k_fd**2 + g2k_fd**2) + phi = cmath.phase(complex(g1k_fd, g2k_fd)) + # g_abs = 0.7 + g1k_fd = g_abs * np.cos(phi - 2*img_rot) + g2k_fd = g_abs * np.sin(phi - 2*img_rot) + # g1k_fd = g_abs * np.cos(phi - 2*img_rot) + # g2k_fd = -g_abs * np.sin(phi - 2*img_rot) + # g1k_fd = g_abs * np.cos(0. - 2*img_rot) + # g2k_fd = -g_abs * np.sin(0. - 2*img_rot) + fd_shear = galsim.Shear(g1=g1k_fd, g2=g2k_fd) return galsim.PositionD(x, y), fd_shear diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index fcdd6f8d3b47edfc9500f455092b4c129c9e3f4f..7121c9e5c3ae851668f5570abca347cf6469fa17 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -341,6 +341,8 @@ class PSFInterp(PSFModel): if findNeighMode == 'hoclistFind': assert(self.hoc != 0), 'hoclist should be built correctly!' imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True) + + imPSf = np.transpose(imPSF) ############TEST: START TestGaussian = False diff --git a/config/config_C6.yaml b/config/config_C6.yaml index ff0e202ea0b84d3fcad52d9246197365da973569..c7b089d4c34f8fbe40c92fbbfc07fc3bd19deffd 100644 --- a/config/config_C6.yaml +++ b/config/config_C6.yaml @@ -9,13 +9,13 @@ # Base diretories and naming setup # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent -work_dir: "/share/home/fangyuedong/fgs_sim/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/sim_v2/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "C6_spectroscopic_20230225" +run_name: "C6_test_profiler" # Whether to use MPI run_option: - use_mpi: YES + use_mpi: NO # NOTE: "n_threads" paramters is currently not used in the backend # simulation codes. It should be implemented later in the web frontend # in order to config the number of threads to request from NAOC cluster @@ -100,7 +100,7 @@ obs_setting: # - give a list of indexes of chips: [ip_1, ip_2...] # - run all chips: null # Note: for all pointings - run_chips: [24] + run_chips: [8] # Whether to enable astrometric modeling enable_astrometric_model: True diff --git a/config/config_C6_test_wcs.yaml b/config/config_C6_test_wcs.yaml new file mode 100644 index 0000000000000000000000000000000000000000..3733c7d72ac3d3f15304695b424df24877246549 --- /dev/null +++ b/config/config_C6_test_wcs.yaml @@ -0,0 +1,222 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# Base diretories and naming setup +# Can add some of the command-line arguments here as well; +# OK to pass either way or both, as long as they are consistent +work_dir: "/share/home/fangyuedong/sim_v2/csst-simulation/workplace/" +data_dir: "/share/simudata/CSSOSDataProductsSims/data/" +run_name: "wcs_test" + +# Whether to use MPI +run_option: + use_mpi: NO + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "Catalog_C6_20221212" + star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + galaxy_cat: "cat2CSSTSim_bundle/" + AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "quickspeclib_ross13.fits" + AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: YES + + # 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 + survey_type: "Photometric" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" + pointing_file: "pointing_radec_246.5_40.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [0] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [9] + + # Whether to enable astrometric modeling + enable_astrometric_model: True + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + # "extra": from seprate file + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + + # Extra shear catalog + # (currently not used) + # shear_cat: "mockShear.cat" + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: OFF # Whether to add field distortions + add_back: OFF # Whether to add sky background + add_dark: OFF # Whether to add dark noise + add_readout: OFF # Whether to add read-out (Gaussian) noise + add_bias: OFF # Whether to add bias-level to images + bias_16channel: OFF # Whether to add different biases for 16 channels + gain_16channel: OFF # Whether to make different gains for 16 channels + shutter_effect: OFF # Whether to add shutter effect + flat_fielding: OFF # Whether to add flat-fielding effect + prnu_effect: OFF # Whether to add PRNU effect + non_linear: OFF # Whether to add non-linearity + cosmic_ray: OFF # Whether to add cosmic-ray + cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: OFF # Whether to simulate CTE trails + saturbloom: OFF # Whether to simulate Saturation & Blooming + add_badcolumns: OFF # Whether to add bad columns + add_hotpixels: OFF # Whether to add hot pixels + add_deadpixels: OFF # Whether to add dead(dark) pixels + bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file diff --git a/profile_C6.sh b/profile_C6.sh new file mode 100755 index 0000000000000000000000000000000000000000..495a778698f05e6f5b9f0bf7529b3591bd32d376 --- /dev/null +++ b/profile_C6.sh @@ -0,0 +1,12 @@ +#!/bin/bash + +date + +python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/sim_v2/csst-simulation/run_sim.py \ + --config_file config_C6_test_wcs.yaml \ + --catalog wcs_test_C6 \ + -c /share/home/fangyuedong/sim_v2/csst-simulation/config + # --config_file test_fd_C6.yaml \ + # --catalog fd_test_C6 \ + # --config_file config_C6.yaml \ + # --catalog C6_Catalog \ diff --git a/run_C6.pbs b/run_C6.pbs index 9fd867219deb03f388f44b297b5bd130f5b8d5c5..dca7a274b69988bb12442c05b43e9e82972b0966 100755 --- a/run_C6.pbs +++ b/run_C6.pbs @@ -1,18 +1,18 @@ #!/bin/bash #PBS -N SIMS -#PBS -l nodes=wcl-1:ppn=60 +#PBS -l nodes=wcl-2:ppn=40 ###PBS -l nodes=wcl-1:ppn=24+wcl-2:ppn=24+wcl-3:ppn=24+wcl-4:ppn=24+wcl-5:ppn=24+wcl-6:ppn=24 #PBS -u fangyuedong ###PBS -j oe cd $PBS_O_WORKDIR -NP=60 +NP=40 date -mpirun -np $NP python3 /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ +mpirun -np $NP python3 /share/home/fangyuedong/sim_v2/csst-simulation/run_sim.py \ --config_file config_C6.yaml \ --catalog C6_Catalog \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config + -c /share/home/fangyuedong/sim_v2/csst-simulation/config diff --git a/test_C6.sh b/test_C6.sh index 82b04e595b71553c1c138f9a1eb5194347677f0e..d976f42dba8cd3185ffeaa61062869a1fc90ee17 100755 --- a/test_C6.sh +++ b/test_C6.sh @@ -2,8 +2,8 @@ date -python3 /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ +python3 /share/home/fangyuedong/sim_v2/csst-simulation/run_sim.py \ --config_file config_C6.yaml \ --catalog C6_Catalog \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config + -c /share/home/fangyuedong/sim_v2/csst-simulation/config