From 96abee7b21defe3bd54a98556812d36fe8ecb9c7 Mon Sep 17 00:00:00 2001
From: fangyuedong <fangyuedong@qxslab.cn>
Date: Mon, 19 Jun 2023 03:23:39 +0800
Subject: [PATCH] wcs definition: x:E, y:-N, add transpose to psf image
---
Catalog/wcs_test_C6.py | 515 ++++++++++++++++++++
ObservationSim/Config/Header/ImageHeader.py | 3 +-
ObservationSim/Instrument/Chip/Chip.py | 2 +-
ObservationSim/Instrument/FocalPlane.py | 13 +-
ObservationSim/MockObject/Galaxy.py | 106 ++--
ObservationSim/MockObject/MockObject.py | 20 +-
ObservationSim/ObservationSim.py | 10 +-
ObservationSim/PSF/FieldDistortion.py | 33 +-
ObservationSim/PSF/PSFInterp.py | 2 +
config/config_C6.yaml | 8 +-
config/config_C6_test_wcs.yaml | 222 +++++++++
profile_C6.sh | 12 +
run_C6.pbs | 8 +-
test_C6.sh | 4 +-
14 files changed, 891 insertions(+), 67 deletions(-)
create mode 100644 Catalog/wcs_test_C6.py
create mode 100644 config/config_C6_test_wcs.yaml
create mode 100755 profile_C6.sh
diff --git a/Catalog/wcs_test_C6.py b/Catalog/wcs_test_C6.py
new file mode 100644
index 0000000..9300e38
--- /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 8a9eafe..8def8ab 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 32b9557..f490506 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 55ef478..ddca78f 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 3db13dc..7874d10 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 0b5ff7f..102191e 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 3c9af78..785370b 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 95dc40e..f7ee0f8 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 fcdd6f8..7121c9e 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 ff0e202..c7b089d 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 0000000..3733c7d
--- /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 0000000..495a778
--- /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 9fd8672..dca7a27 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 82b04e5..d976f42 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
--
GitLab