Commit f5c20f83 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge remote-tracking branch 'origin/sim_scheduler' into develop

parents 3a9baf79 53592d2e
...@@ -6,3 +6,5 @@ dist/* ...@@ -6,3 +6,5 @@ dist/*
*disperse.c *disperse.c
*interp.c *interp.c
!*libshao.so !*libshao.so
*.out
pnodes
\ No newline at end of file
...@@ -39,9 +39,9 @@ bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5', ...@@ -39,9 +39,9 @@ bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5',
qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits'] qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits']
# star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5'] # star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5']
star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)] star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.),[246.5, 40]]
star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5'] star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5','trilegal_calMag_mpi_Nside_128_healpix.hdf5']
class StarParm(ctypes.Structure): class StarParm(ctypes.Structure):
_fields_ = [ _fields_ = [
...@@ -187,7 +187,7 @@ class Catalog(CatalogBase): ...@@ -187,7 +187,7 @@ class Catalog(CatalogBase):
rv_c = obj.param['rv']/(atcons.c.value/1000.) rv_c = obj.param['rv']/(atcons.c.value/1000.)
Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c)) Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c))
wave_RV = wave*Doppler_factor wave_RV = wave*Doppler_factor
return wave_RV, spec return wave_RV, np.power(10,spec[:])
def _load_SED_lib_gals(self): def _load_SED_lib_gals(self):
pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
...@@ -511,7 +511,6 @@ class Catalog(CatalogBase): ...@@ -511,7 +511,6 @@ class Catalog(CatalogBase):
# erg/s/cm2/A --> photon/s/m2/A # erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
if obj.type == 'quasar': if obj.type == 'quasar':
# integrate to get the magnitudes # integrate to get the magnitudes
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
......
...@@ -3,8 +3,9 @@ import logging ...@@ -3,8 +3,9 @@ import logging
import ObservationSim.Config._util as _util import ObservationSim.Config._util as _util
from ObservationSim.Config.Header import generatePrimaryHeader from ObservationSim.Config.Header import generatePrimaryHeader
class ChipOutput(object): class ChipOutput(object):
def __init__(self, config, chip, filt, pointing): def __init__(self, config, chip, filt, pointing, logger_filename=None):
self.config = config self.config = config
self.chip = chip self.chip = chip
self.filt = filt self.filt = filt
...@@ -15,12 +16,12 @@ class ChipOutput(object): ...@@ -15,12 +16,12 @@ class ChipOutput(object):
self.h_prim = generatePrimaryHeader( self.h_prim = generatePrimaryHeader(
xlen=chip.npix_x, xlen=chip.npix_x,
ylen=chip.npix_y, ylen=chip.npix_y,
pointing_id = pointing.obs_id, pointing_id=pointing.obs_id,
pointing_type_code = pointing.pointing_type_code, pointing_type_code=pointing.pointing_type_code,
ra=pointing.ra, ra=pointing.ra,
dec=pointing.dec, dec=pointing.dec,
pixel_scale=chip.pix_scale, pixel_scale=chip.pix_scale,
time_pt = pointing.timestamp, time_pt=pointing.timestamp,
exptime=pointing.exp_time, exptime=pointing.exp_time,
im_type=pointing.pointing_type, im_type=pointing.pointing_type,
sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
...@@ -29,18 +30,22 @@ class ChipOutput(object): ...@@ -29,18 +30,22 @@ class ChipOutput(object):
run_counter=self.config["run_counter"], run_counter=self.config["run_counter"],
chip_name=self.chip_label) chip_name=self.chip_label)
obs_id = _util.get_obs_id(img_type=self.pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_id=pointing.obs_id, pointing_type_code = pointing.pointing_type_code) obs_id = _util.get_obs_id(img_type=self.pointing_type, project_cycle=config["project_cycle"], run_counter=config[
"run_counter"], pointing_id=pointing.obs_id, pointing_type_code=pointing.pointing_type_code)
self.subdir = pointing.output_dir self.subdir = pointing.output_dir
self.cat_name = self.h_prim['FILENAME'] + '.cat' self.cat_name = self.h_prim['FILENAME'] + '.cat'
if logger_filename is None:
logger_filename = self.h_prim['FILENAME'] + '.log' logger_filename = self.h_prim['FILENAME'] + '.log'
self.logger = logging.getLogger() self.logger = logging.getLogger()
fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8') fh = logging.FileHandler(os.path.join(
self.subdir, logger_filename), mode='w+', encoding='utf-8')
fh.setLevel(logging.DEBUG) fh.setLevel(logging.DEBUG)
self.logger.setLevel(logging.DEBUG) self.logger.setLevel(logging.DEBUG)
logging.getLogger('numba').setLevel(logging.WARNING) logging.getLogger('numba').setLevel(logging.WARNING)
formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s') formatter = logging.Formatter(
'%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s')
fh.setFormatter(formatter) fh.setFormatter(formatter)
self.logger.addHandler(fh) self.logger.addHandler(fh)
...@@ -51,7 +56,7 @@ class ChipOutput(object): ...@@ -51,7 +56,7 @@ class ChipOutput(object):
self.hdr = hdr1 + hdr2 self.hdr = hdr1 + hdr2
self.fmt = fmt1 + fmt2 self.fmt = fmt1 + fmt2
self.logger.info("pointing_type = %s\n"%(self.pointing_type)) self.logger.info("pointing_type = %s\n" % (self.pointing_type))
def Log_info(self, message): def Log_info(self, message):
print(message) print(message)
...@@ -67,7 +72,8 @@ class ChipOutput(object): ...@@ -67,7 +72,8 @@ class ChipOutput(object):
def create_output_file(self): def create_output_file(self):
if self.pointing_type == 'SCI': if self.pointing_type == 'SCI':
self.cat = open(os.path.join(self.subdir, self.cat_name), "w") self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) self.logger.info("Creating catalog file %s ...\n" %
(os.path.join(self.subdir, self.cat_name)))
if not self.hdr.endswith("\n"): if not self.hdr.endswith("\n"):
self.hdr += "\n" self.hdr += "\n"
self.cat.write(self.hdr) self.cat.write(self.hdr)
...@@ -76,8 +82,9 @@ class ChipOutput(object): ...@@ -76,8 +82,9 @@ class ChipOutput(object):
ximg = obj.real_pos.x + 1.0 ximg = obj.real_pos.x + 1.0
yimg = obj.real_pos.y + 1.0 yimg = obj.real_pos.y + 1.0
line = self.fmt%( line = self.fmt % (
obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(
self.filt), obj.type,
obj.pmra, obj.pmdec, obj.rv, obj.parallax) obj.pmra, obj.pmdec, obj.rv, obj.parallax)
line += obj.additional_output_str line += obj.additional_output_str
if not line.endswith("\n"): if not line.endswith("\n"):
......
...@@ -27,52 +27,52 @@ class CatalogBase(metaclass=ABCMeta): ...@@ -27,52 +27,52 @@ class CatalogBase(metaclass=ABCMeta):
@staticmethod @staticmethod
def initialize_param(): def initialize_param():
param = { param = {
"star":-1, "star": -1,
"id":-1, "id": -1,
"ra":0, "ra": 0,
"dec":0., "dec": 0.,
"ra_orig":0., "ra_orig": 0.,
"dec_orig":0., "dec_orig": 0.,
"z":0., "z": 0.,
"sed_type":-1, "sed_type": -1,
"model_tag":"unknown", "model_tag": "unknown",
"mag_use_normal":100., "mag_use_normal": 100.,
"theta":0., "theta": 0.,
"kappa":0., "kappa": 0.,
"g1":0., "g1": 0.,
"g2":0., "g2": 0.,
"bfrac":0., "bfrac": 0.,
"av":0., "av": 0.,
"redden":0., "redden": 0.,
"hlr_bulge":0., "hlr_bulge": 0.,
"hlr_disk":0., "hlr_disk": 0.,
"ell_bulge":0., "ell_bulge": 0.,
"ell_disk":0., "ell_disk": 0.,
"ell_tot":0., "ell_tot": 0.,
"e1_disk":0., "e1_disk": 0.,
"e2_disk":0., "e2_disk": 0.,
"e1_bulge":0., "e1_bulge": 0.,
"e2_bulge":0., "e2_bulge": 0.,
"teff":0., "teff": 0.,
"logg":0., "logg": 0.,
"feh":0., "feh": 0.,
"DM":0., "DM": 0.,
"stellarMass":1., "stellarMass": 1.,
# C6 galaxies parameters # C6 galaxies parameters
"e1":0., "e1": 0.,
"e2":0., "e2": 0.,
"bulgemass":0., "bulgemass": 0.,
"diskmass":0., "diskmass": 0.,
"size":0., "size": 0.,
"detA":0., "detA": 0.,
"type":0, "type": 0,
"veldisp":0., "veldisp": 0.,
"coeff": np.zeros(20), "coeff": np.zeros(20),
# Astrometry related # Astrometry related
"pmra":0., "pmra": 0.,
"pmdec":0., "pmdec": 0.,
"rv":0., "rv": 0.,
"parallax":1e-9 "parallax": 1e-9
} }
return param return param
...@@ -87,27 +87,32 @@ class CatalogBase(metaclass=ABCMeta): ...@@ -87,27 +87,32 @@ class CatalogBase(metaclass=ABCMeta):
return e1, e2, e_total return e1, e2, e_total
@staticmethod @staticmethod
def convert_sed(mag, sed, target_filt, norm_filt=None): def convert_sed(mag, sed, target_filt, norm_filt=None, mu=1.):
bandpass = target_filt.bandpass_full bandpass = target_filt.bandpass_full
if norm_filt is not None: if norm_filt is not None:
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001 norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
else: else:
norm_filt = Table( norm_filt = Table(
np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY']) np.array(np.array([bandpass.wave_list*10.0, bandpass.func(
bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
) )
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001 norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag, sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
spectrum=sed, spectrum=sed,
norm_thr=norm_filt, norm_thr=norm_filt,
sWave=np.floor(norm_filt[norm_thr_rang_ids][0][0]), sWave=np.floor(
norm_filt[norm_thr_rang_ids][0][0]),
eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0])) eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
sed_photon = copy.copy(sed) sed_photon = copy.copy(sed)
sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T sed_photon = np.array(
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') [sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1] * mu), interpolant='nearest')
# Get magnitude # Get magnitude
sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) sed_photon = galsim.SED(sed_photon, wave_type='A',
flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass) interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
mag_csst = getABMAG( mag_csst = getABMAG(
interFlux=interFlux, interFlux=interFlux,
......
...@@ -8,6 +8,7 @@ from ObservationSim.MockObject.MockObject import MockObject ...@@ -8,6 +8,7 @@ from ObservationSim.MockObject.MockObject import MockObject
# import tracemalloc # import tracemalloc
class Galaxy(MockObject): class Galaxy(MockObject):
def __init__(self, param, logger=None): def __init__(self, param, logger=None):
super().__init__(param, logger=logger) super().__init__(param, logger=logger)
...@@ -16,6 +17,11 @@ class Galaxy(MockObject): ...@@ -16,6 +17,11 @@ class Galaxy(MockObject):
self.disk_sersic_idx = 1. self.disk_sersic_idx = 1.
if not hasattr(self, "bulge_sersic_idx"): if not hasattr(self, "bulge_sersic_idx"):
self.bulge_sersic_idx = 4. self.bulge_sersic_idx = 4.
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def unload_SED(self): def unload_SED(self):
"""(Test) free up SED memory """(Test) free up SED memory
...@@ -24,14 +30,16 @@ class Galaxy(MockObject): ...@@ -24,14 +30,16 @@ class Galaxy(MockObject):
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None): def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
if len(psf_list) != len(bandpass_list): if len(psf_list) != len(bandpass_list):
raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.") raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = [] objs = []
if nphotons_tot == None: if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime) nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot) # print("nphotons_tot = ", nphotons_tot)
try: try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e: except Exception as e:
print(e) print(e)
if self.logger: if self.logger:
...@@ -54,10 +62,12 @@ class Galaxy(MockObject): ...@@ -54,10 +62,12 @@ class Galaxy(MockObject):
return -1 return -1
psf = psf_list[i] psf = psf_list[i]
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0) disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape) disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0) bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape) bulge = bulge.shear(bulge_shape)
...@@ -67,13 +77,15 @@ class Galaxy(MockObject): ...@@ -67,13 +77,15 @@ class Galaxy(MockObject):
gal = bulge gal = bulge
else: else:
gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal = gal.withFlux(nphotons)
if fd_shear is not None: if fd_shear is not None:
g1 += fd_shear.g1 g1 += fd_shear.g1
g2 += fd_shear.g2 g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2) gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear) gal = gal.shear(gal_shear)
# Magnification
gal = gal.magnify(self.mu)
gal = galsim.Convolve(psf, gal) gal = galsim.Convolve(psf, gal)
gal = gal.withFlux(nphotons)
objs.append(gal) objs.append(gal)
final = galsim.Sum(objs) final = galsim.Sum(objs)
...@@ -85,7 +97,8 @@ class Galaxy(MockObject): ...@@ -85,7 +97,8 @@ class Galaxy(MockObject):
# print("nphotons_tot = ", nphotons_tot) # print("nphotons_tot = ", nphotons_tot)
try: try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e: except Exception as e:
print(e) print(e)
if self.logger: if self.logger:
...@@ -121,10 +134,12 @@ class Galaxy(MockObject): ...@@ -121,10 +134,12 @@ class Galaxy(MockObject):
is_updated = 0 is_updated = 0
# Model the galaxy as disk + bulge # Model the galaxy as disk + bulge
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp) 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_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape) 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 = 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_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape) bulge = bulge.shear(bulge_shape)
...@@ -155,7 +170,8 @@ class Galaxy(MockObject): ...@@ -155,7 +170,8 @@ class Galaxy(MockObject):
# print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
# Get PSF model # Get PSF model
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) psf, pos_shear = psf_model.get_PSF(
chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
if self.bfrac == 0: if self.bfrac == 0:
gal_temp = disk gal_temp = disk
...@@ -164,10 +180,13 @@ class Galaxy(MockObject): ...@@ -164,10 +180,13 @@ class Galaxy(MockObject):
else: else:
gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal_temp = gal_temp.shear(gal_shear) gal_temp = gal_temp.shear(gal_shear)
gal_temp = gal_temp.withFlux(nphotons) # Magnification
gal_temp = gal_temp.magnify(self.mu)
if not big_galaxy: # Not apply PSF for very big galaxy if not big_galaxy: # Not apply PSF for very big galaxy
gal_temp = galsim.Convolve(psf, gal_temp) gal_temp = galsim.Convolve(psf, gal_temp)
gal_temp = gal_temp.withFlux(nphotons)
if i == 0: if i == 0:
gal = gal_temp gal = gal_temp
else: else:
...@@ -184,7 +203,8 @@ class Galaxy(MockObject): ...@@ -184,7 +203,8 @@ class Galaxy(MockObject):
# ERROR happens # ERROR happens
return 2, pos_shear return 2, pos_shear
stamp.setCenter(x_nominal, y_nominal) stamp.setCenter(x_nominal, y_nominal)
bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) bounds = stamp.bounds & galsim.BoundsI(
0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() > 0: if bounds.area() > 0:
chip.img.setOrigin(0, 0) chip.img.setOrigin(0, 0)
chip.img[bounds] += stamp[bounds] chip.img[bounds] += stamp[bounds]
...@@ -194,9 +214,9 @@ class Galaxy(MockObject): ...@@ -194,9 +214,9 @@ class Galaxy(MockObject):
if is_updated == 0: if is_updated == 0:
# Return code 0: object photons missed this detector # Return code 0: object photons missed this detector
print("obj %s missed"%(self.id)) print("obj %s missed" % (self.id))
if self.logger: if self.logger:
self.logger.info("obj %s missed"%(self.id)) self.logger.info("obj %s missed" % (self.id))
return 0, pos_shear return 0, pos_shear
# # [C6 TEST] # # [C6 TEST]
...@@ -209,7 +229,8 @@ class Galaxy(MockObject): ...@@ -209,7 +229,8 @@ class Galaxy(MockObject):
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed, sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
norm_thr=normFilter, norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]), sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
if sedNormFactor == 0: if sedNormFactor == 0:
return 2, None return 2, None
...@@ -230,7 +251,6 @@ class Galaxy(MockObject): ...@@ -230,7 +251,6 @@ class Galaxy(MockObject):
chip_wcs_local = self.chip_wcs.local(self.real_pos) chip_wcs_local = self.chip_wcs.local(self.real_pos)
big_galaxy = False big_galaxy = False
if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
big_galaxy = True big_galaxy = True
...@@ -244,7 +264,8 @@ class Galaxy(MockObject): ...@@ -244,7 +264,8 @@ class Galaxy(MockObject):
flat_cube = chip.flat_cube flat_cube = chip.flat_cube
xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062,
'C': 4.035447379743442, 'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos grating_split_pos_chip = 0 + grating_split_pos
branges = np.zeros([len(bandpass_list), 2]) branges = np.zeros([len(bandpass_list), 2])
...@@ -267,10 +288,12 @@ class Galaxy(MockObject): ...@@ -267,10 +288,12 @@ class Galaxy(MockObject):
brange = branges[i] brange = branges[i]
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) # 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 = 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_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape) 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 = 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_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape) bulge = bulge.shear(bulge_shape)
...@@ -286,12 +309,13 @@ class Galaxy(MockObject): ...@@ -286,12 +309,13 @@ class Galaxy(MockObject):
# kfrac = np.random.random()*(1.0 - self.bfrac) # kfrac = np.random.random()*(1.0 - self.bfrac)
# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
gal = gal.withFlux(tel.pupil_area * exptime)
if fd_shear: if fd_shear:
g1 += fd_shear.g1 g1 += fd_shear.g1
g2 += fd_shear.g2 g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2) gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear) gal = gal.shear(gal_shear)
gal = gal.magnify(self.mu)
gal = gal.withFlux(tel.pupil_area * exptime)
# gal = galsim.Convolve(psf, gal) # gal = galsim.Convolve(psf, gal)
# if not big_galaxy: # Not apply PSF for very big galaxy # if not big_galaxy: # Not apply PSF for very big galaxy
...@@ -299,23 +323,25 @@ class Galaxy(MockObject): ...@@ -299,23 +323,25 @@ class Galaxy(MockObject):
# # if fd_shear is not None: # # if fd_shear is not None:
# # gal = gal.shear(fd_shear) # # gal = gal.shear(fd_shear)
starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space') starImg = gal.drawImage(
wcs=chip_wcs_local, offset=offset, method='real_space')
origin_star = [y_nominal - (starImg.center.y - starImg.ymin), origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)] x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0, 0) starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]] gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1] gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]: if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse # part img disperse
subImg_p1 = starImg.array[:, 0:subSlitPos] subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1) star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0) star_p1.setOrigin(0, 0)
origin_p1 = origin_star origin_p1 = origin_star
xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0 xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-0 ycenter_p1 = y_nominal-0
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
...@@ -330,9 +356,10 @@ class Galaxy(MockObject): ...@@ -330,9 +356,10 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] subImg_p2 = starImg.array[:,
subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0) star_p2.setOrigin(0, 0)
origin_p2 = [origin_star[0], grating_split_pos_chip] origin_p2 = [origin_star[0], grating_split_pos_chip]
...@@ -351,11 +378,11 @@ class Galaxy(MockObject): ...@@ -351,11 +378,11 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1 del sdp_p1
del sdp_p2 del sdp_p2
elif grating_split_pos_chip<=gal_origin[1]: elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
...@@ -367,9 +394,9 @@ class Galaxy(MockObject): ...@@ -367,9 +394,9 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
elif grating_split_pos_chip>=gal_end[1]: elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
...@@ -381,7 +408,7 @@ class Galaxy(MockObject): ...@@ -381,7 +408,7 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin) # print(self.y_nominal, starImg.center.y, starImg.ymin)
...@@ -391,11 +418,13 @@ class Galaxy(MockObject): ...@@ -391,11 +418,13 @@ class Galaxy(MockObject):
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.): def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
if flux == None: if flux == None:
flux = self.getElectronFluxFilt(filt, tel, exptime) flux = self.getElectronFluxFilt(filt, tel, exptime)
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0) disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk) disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape) disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0) bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape) bulge = bulge.shear(bulge_shape)
...@@ -407,5 +436,6 @@ class Galaxy(MockObject): ...@@ -407,5 +436,6 @@ class Galaxy(MockObject):
return final return final
def getObservedEll(self, g1=0, g2=0): def getObservedEll(self, g1=0, g2=0):
e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2) e1_obs, e2_obs, e_obs, theta = eObs(
self.e1_total, self.e2_total, g1, g2)
return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs
import galsim import galsim
import os, sys import os
import sys
import numpy as np import numpy as np
import astropy.constants as cons import astropy.constants as cons
from astropy.table import Table from astropy.table import Table
...@@ -8,9 +9,15 @@ from scipy import interpolate ...@@ -8,9 +9,15 @@ from scipy import interpolate
from ObservationSim.MockObject.MockObject import MockObject from ObservationSim.MockObject.MockObject import MockObject
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
class Quasar(MockObject): class Quasar(MockObject):
def __init__(self, param, logger=None): def __init__(self, param, logger=None):
super().__init__(param, logger=logger) super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
''' '''
...@@ -22,8 +29,9 @@ class Quasar(MockObject): ...@@ -22,8 +29,9 @@ class Quasar(MockObject):
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
if sed_templates is None: if sed_templates is None:
# Read SED data directly # Read SED data directly
itype = objtypes[cosids==self.sed_type][0] itype = objtypes[cosids == self.sed_type][0]
sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type)) sed_file = os.path.join(
sed_path, itype + "_ID%s.sed" % (self.sed_type))
if not os.path.exists(sed_file): if not os.path.exists(sed_file):
raise ValueError("!!! No SED found.") raise ValueError("!!! No SED found.")
sed_data = Table.read(sed_file, format="ascii") sed_data = Table.read(sed_file, format="ascii")
...@@ -40,20 +48,26 @@ class Quasar(MockObject): ...@@ -40,20 +48,26 @@ class Quasar(MockObject):
wave, flux = sed_data[0], sed_data[1] wave, flux = sed_data[0], sed_data[1]
flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13 flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX')) sed_photon = Table(
np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
# Get scaling factor for SED # Get scaling factor for SED
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
spectrum=sed_photon, spectrum=sed_photon,
norm_thr=normFilter, norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]), sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T sed_photon = np.array(
[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
# Convert to galsim.SED object # Convert to galsim.SED object
spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False) sed_photon[:, 1]), interpolant='nearest')
self.sed = galsim.SED(spec, wave_type='A',
flux_type='1', fast=False)
# Get magnitude # Get magnitude
interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full) interFlux = integrate_sed_bandpass(
self.param['mag_%s'%target_filt.filter_type] = getABMAG( sed=self.sed, bandpass=target_filt.bandpass_full)
self.param['mag_%s' % target_filt.filter_type] = getABMAG(
interFlux=interFlux, interFlux=interFlux,
bandpass=target_filt.bandpass_full) bandpass=target_filt.bandpass_full)
# print('mag_use_normal = ', self.param['mag_use_normal']) # print('mag_use_normal = ', self.param['mag_use_normal'])
...@@ -61,7 +75,8 @@ class Quasar(MockObject): ...@@ -61,7 +75,8 @@ class Quasar(MockObject):
elif survey_type == "spectroscopic": elif survey_type == "spectroscopic":
if sed_templates is None: if sed_templates is None:
self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes) self.sedPhotons(sed_path=sed_path,
cosids=cosids, objtypes=objtypes)
else: else:
sed_data = sed_templates[self.sed_type] sed_data = sed_templates[self.sed_type]
sed_data = getObservedSED( sed_data = getObservedSED(
...@@ -74,8 +89,8 @@ class Quasar(MockObject): ...@@ -74,8 +89,8 @@ class Quasar(MockObject):
y = speci(lamb) y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A # erg/s/cm2/A --> photo/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) self.sed = Table(
np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
def unload_SED(self): def unload_SED(self):
"""(Test) free up SED memory """(Test) free up SED memory
......
import galsim import galsim
import os, sys import os
import sys
import numpy as np import numpy as np
import astropy.constants as cons import astropy.constants as cons
from astropy.table import Table from astropy.table import Table
...@@ -8,9 +9,12 @@ from scipy import interpolate ...@@ -8,9 +9,12 @@ from scipy import interpolate
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed
from ObservationSim.MockObject.MockObject import MockObject from ObservationSim.MockObject.MockObject import MockObject
class Star(MockObject): class Star(MockObject):
def __init__(self, param, logger=None): def __init__(self, param, logger=None):
super().__init__(param, logger=logger) super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
self.mu = 1.
def unload_SED(self): def unload_SED(self):
"""(Test) free up SED memory """(Test) free up SED memory
...@@ -28,13 +32,15 @@ class Star(MockObject): ...@@ -28,13 +32,15 @@ class Star(MockObject):
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.): def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
if len(psf_list) != len(bandpass_list): if len(psf_list) != len(bandpass_list):
raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.") raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = [] objs = []
if nphotons_tot == None: if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime) nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
try: try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full) full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e: except Exception as e:
print(e) print(e)
self.logger.error(e) self.logger.error(e)
......
...@@ -10,6 +10,7 @@ from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSL ...@@ -10,6 +10,7 @@ from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSL
from astropy.time import Time from astropy.time import Time
from datetime import datetime, timezone from datetime import datetime, timezone
def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get exposure time # Get exposure time
...@@ -20,9 +21,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -20,9 +21,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Load catalogues # Load catalogues
if catalog is None: if catalog is None:
self.chip_output.Log_error("Catalog interface class must be specified for SCIE-OBS") self.chip_output.Log_error(
raise ValueError("Catalog interface class must be specified for SCIE-OBS") "Catalog interface class must be specified for SCIE-OBS")
cat = catalog(config=self.overall_config, chip=chip, pointing=pointing, chip_output=self.chip_output, filt=filt) raise ValueError(
"Catalog interface class must be specified for SCIE-OBS")
cat = catalog(config=self.overall_config, chip=chip,
pointing=pointing, chip_output=self.chip_output, filt=filt)
# Prepare output file(s) for this chip # Prepare output file(s) for this chip
# [NOTE] Headers of output .cat file may be updated by Catalog instance # [NOTE] Headers of output .cat file may be updated by Catalog instance
...@@ -31,12 +35,15 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -31,12 +35,15 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Prepare the PSF model # Prepare the PSF model
if self.overall_config["psf_setting"]["psf_model"] == "Gauss": if self.overall_config["psf_setting"]["psf_model"] == "Gauss":
psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"]) psf_model = PSFGauss(
chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"])
elif self.overall_config["psf_setting"]["psf_model"] == "Interp": elif self.overall_config["psf_setting"]["psf_model"] == "Interp":
if chip.survey_type == "spectroscopic": if chip.survey_type == "spectroscopic":
psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"]) psf_model = PSFInterpSLS(
chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"])
else: else:
psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"]) psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples,
PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"])
else: else:
self.chip_output.Log_error("unrecognized PSF model type!!", flush=True) self.chip_output.Log_error("unrecognized PSF model type!!", flush=True)
...@@ -50,20 +57,22 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -50,20 +57,22 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get the filter which will be used for magnitude cut # Get the filter which will be used for magnitude cut
for ifilt in range(len(self.all_filters)): for ifilt in range(len(self.all_filters)):
temp_filter = self.all_filters[ifilt] temp_filter = self.all_filters[ifilt]
temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) temp_filter.update_limit_saturation_mags(
exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip)
if temp_filter.filter_type.lower() == self.overall_config["obs_setting"]["cut_in_band"].lower(): if temp_filter.filter_type.lower() == self.overall_config["obs_setting"]["cut_in_band"].lower():
cut_filter = temp_filter cut_filter = temp_filter
# Read in shear values from configuration file if the constant shear type is used # Read in shear values from configuration file if the constant shear type is used
if self.overall_config["shear_setting"]["shear_type"] == "constant": if self.overall_config["shear_setting"]["shear_type"] == "constant":
g1_field, g2_field, _ = get_shear_field(config=self.overall_config) g1_field, g2_field, _ = get_shear_field(config=self.overall_config)
self.chip_output.Log_info("Use constant shear: g1=%.5f, g2=%.5f"%(g1_field, g2_field)) self.chip_output.Log_info(
"Use constant shear: g1=%.5f, g2=%.5f" % (g1_field, g2_field))
# Get chip WCS # Get chip WCS
if not hasattr(self, 'h_ext'): if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing) _, _ = self.prepare_headers(chip=chip, pointing=pointing)
chip_wcs = galsim.FitsWCS(header = self.h_ext) chip_wcs = galsim.FitsWCS(header=self.h_ext)
# Loop over objects # Loop over objects
nobj = len(cat.objs) nobj = len(cat.objs)
...@@ -80,17 +89,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -80,17 +89,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
try: try:
sed_data = cat.load_sed(obj) sed_data = cat.load_sed(obj)
norm_filt = cat.load_norm_filt(obj) norm_filt = cat.load_norm_filt(obj)
obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = cat.convert_sed( obj.sed, obj.param["mag_%s" % filt.filter_type.lower()], obj.param["flux_%s" % filt.filter_type.lower()] = cat.convert_sed(
mag=obj.param["mag_use_normal"], mag=obj.param["mag_use_normal"],
sed=sed_data, sed=sed_data,
target_filt=filt, target_filt=filt,
norm_filt=norm_filt, norm_filt=norm_filt,
mu=obj.mu
) )
_, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = cat.convert_sed( _, obj.param["mag_%s" % cut_filter.filter_type.lower()], obj.param["flux_%s" % cut_filter.filter_type.lower()] = cat.convert_sed(
mag=obj.param["mag_use_normal"], mag=obj.param["mag_use_normal"],
sed=sed_data, sed=sed_data,
target_filt=cut_filter, target_filt=cut_filter,
norm_filt=norm_filt, norm_filt=norm_filt,
mu=obj.mu
) )
except Exception as e: except Exception as e:
traceback.print_exc() traceback.print_exc()
...@@ -102,16 +113,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -102,16 +113,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Exclude very bright/dim objects (for now) # Exclude very bright/dim objects (for now)
if cut_filter.is_too_bright( if cut_filter.is_too_bright(
mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()], mag=obj.param["mag_%s" %
self.overall_config["obs_setting"]["cut_in_band"].lower()],
margin=self.overall_config["obs_setting"]["mag_sat_margin"]): margin=self.overall_config["obs_setting"]["mag_sat_margin"]):
self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()])) self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f" % (
obj.id, cut_filter.filter_type, obj.param["mag_%s" % self.overall_config["obs_setting"]["cut_in_band"].lower()]))
bright_obj += 1 bright_obj += 1
obj.unload_SED() obj.unload_SED()
continue continue
if filt.is_too_dim( if filt.is_too_dim(
mag=obj.getMagFilter(filt), mag=obj.getMagFilter(filt),
margin=self.overall_config["obs_setting"]["mag_lim_margin"]): margin=self.overall_config["obs_setting"]["mag_lim_margin"]):
self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt))) self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f" % (
obj.id, filt.filter_type, obj.getMagFilter(filt)))
dim_obj += 1 dim_obj += 1
obj.unload_SED() obj.unload_SED()
continue continue
...@@ -130,14 +144,16 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -130,14 +144,16 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
raise ValueError("Unknown shear input") raise ValueError("Unknown shear input")
# Get position of object on the focal plane # Get position of object on the focal plane
pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext) pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(
img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext)
# [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # [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 # Otherwise they will be considered missed objects
# if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)): # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)):
if pos_img.x == -1 or pos_img.y == -1: if pos_img.x == -1 or pos_img.y == -1:
self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f' % (
self.chip_output.Log_error("Objected missed: %s"%(obj.id)) obj.ra, obj.dec, obj.ra_orig, obj.dec_orig))
self.chip_output.Log_error("Objected missed: %s" % (obj.id))
missed_obj += 1 missed_obj += 1
obj.unload_SED() obj.unload_SED()
continue continue
...@@ -146,7 +162,8 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -146,7 +162,8 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
try: try:
if self.overall_config["run_option"]["out_cat_only"]: if self.overall_config["run_option"]["out_cat_only"]:
isUpdated = True isUpdated = True
obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) obj.real_pos = obj.getRealPos(
chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs)
pos_shear = 0. pos_shear = 0.
elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]: elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]:
isUpdated, pos_shear = obj.drawObj_multiband( isUpdated, pos_shear = obj.drawObj_multiband(
...@@ -181,9 +198,10 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -181,9 +198,10 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
pass pass
elif isUpdated == 0: elif isUpdated == 0:
missed_obj += 1 missed_obj += 1
self.chip_output.Log_error("Objected missed: %s"%(obj.id)) self.chip_output.Log_error("Objected missed: %s" % (obj.id))
else: else:
self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id)) self.chip_output.Log_error(
"Draw error, object omitted: %s" % (obj.id))
continue continue
except Exception as e: except Exception as e:
traceback.print_exc() traceback.print_exc()
...@@ -196,38 +214,47 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -196,38 +214,47 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
del psf_model del psf_model
gc.collect() gc.collect()
self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB" %
(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)))
self.chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, nobj)) self.chip_output.Log_info(
self.chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, nobj)) "# objects that are too bright %d out of %d" % (bright_obj, nobj))
self.chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, nobj)) self.chip_output.Log_info(
"# objects that are too dim %d out of %d" % (dim_obj, nobj))
self.chip_output.Log_info(
"# objects that are missed %d out of %d" % (missed_obj, nobj))
# Apply flat fielding (with shutter effects) # Apply flat fielding (with shutter effects)
flat_normal = np.ones_like(chip.img.array) flat_normal = np.ones_like(chip.img.array)
if obs_param["flat_fielding"] == True: if obs_param["flat_fielding"] == True:
flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array) flat_normal = flat_normal * chip.flat_img.array / \
np.mean(chip.flat_img.array)
if obs_param["shutter_effect"] == True: if obs_param["shutter_effect"] == True:
flat_normal = flat_normal * chip.shutter_img flat_normal = flat_normal * chip.shutter_img
flat_normal = np.array(flat_normal, dtype='float32') flat_normal = np.array(flat_normal, dtype='float32')
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True]) self.updateHeaderInfo(header_flag='ext', keys=[
'SHTSTAT'], values=[True])
else: else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']]) self.updateHeaderInfo(header_flag='ext', keys=['SHTSTAT', 'SHTOPEN1', 'SHTCLOS0'], values=[
True, self.h_ext['SHTCLOS1'], self.h_ext['SHTOPEN0']])
chip.img *= flat_normal chip.img *= flat_normal
del flat_normal del flat_normal
# renew header info # renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs) t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.s,开始曝光 # ccd刷新2s,等待0.s,开始曝光
t_obs_renew = Time(t_obs.mjd - (2.+0.) / 86400., format="mjd") t_obs_renew = Time(t_obs.mjd - (2.+0.) / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys=[
'DATE-OBS'], values=[t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s # dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.+pointing.exp_time]) self.updateHeaderInfo(header_flag='ext', keys=[
'DARKTIME'], values=[0.+0.+pointing.exp_time])
return chip, filt, tel, pointing return chip, filt, tel, pointing
...@@ -136,7 +136,7 @@ def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -136,7 +136,7 @@ def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time]) self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.0+0.0+pointing.exp_time])
return chip, filt, tel, pointing return chip, filt, tel, pointing
\ No newline at end of file
...@@ -76,7 +76,7 @@ with open("requirements.txt", "r") as f: ...@@ -76,7 +76,7 @@ with open("requirements.txt", "r") as f:
] ]
setup(name='csst_msc_sim', setup(name='csst_msc_sim',
version='3.0.0', version='3.0.0rc',
packages=find_packages(), packages=find_packages(),
# install_requires=[ # install_requires=[
# # 'numpy>=1.18.5', # # 'numpy>=1.18.5',
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment