import unittest from scipy import interpolate import astropy.constants as cons from astropy.table import Table import h5py as h5 import sys import os import math from itertools import islice import numpy as np import galsim import yaml import copy from astropy.cosmology import FlatLambdaCDM from astropy import constants from astropy import units as U from ObservationSim.MockObject._util import getObservedSED from ObservationSim.MockObject import CatalogBase, Galaxy from ObservationSim.Instrument import Chip, Filter, FilterParam, FocalPlane from ObservationSim.PSF.PSFInterp import PSFInterp from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getABMAG class Catalog(CatalogBase): def __init__(self): super().__init__() self.rotation = 0. def load_norm_filt(self, obj): return None def load_sed(self, obj, **kward): pcs = h5.File(os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_msc_sim/csst_fz_msc/sedlibs/pcs.h5'), "r") lamb = h5.File(os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_msc_sim/csst_fz_msc/sedlibs/lamb.h5'), "r") lamb_gal = lamb['lamb'][()] pcs = pcs['pcs'][()] cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) factor = 10**(-.4 * cosmo.distmod(obj.param['z']).value) flux = np.matmul(pcs, obj.param['coeff']) * factor # if np.any(flux < 0): # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) flux[flux < 0] = 0. sedcat = np.vstack((lamb_gal, flux)).T sed_data = getObservedSED( sedCat=sedcat, redshift=obj.param['z'], av=obj.param["av"], redden=obj.param["redden"] ) wave, flux = sed_data[0], sed_data[1] 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')) # obj.param["sed"] = sed del wave del flux return sed def _load(self, file_path): gals_cat = h5.File(file_path, 'r')['galaxies'] for ikeys in gals_cat.keys(): gals = gals_cat[ikeys] param = self.initialize_param() igals = 9 param['z'] = gals['redshift'][igals] param['mag_use_normal'] = gals['mag_csst_g'][igals] print(param['mag_use_normal']) param['e1'] = gals['ellipticity_true'][igals][0] param['e2'] = gals['ellipticity_true'][igals][1] # For shape calculation param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity( e1=gals['ellipticity_true'][igals][0], e2=gals['ellipticity_true'][igals][1], rotation=self.rotation, unit='radians') param['e1_disk'] = param['e1'] param['e2_disk'] = param['e2'] param['e1_bulge'] = param['e1'] param['e2_bulge'] = param['e2'] # Masses param['bulgemass'] = gals['bulgemass'][igals] param['diskmass'] = gals['diskmass'][igals] param['size'] = gals['size'][igals] # 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] # TEST no redening and no extinction param['av'] = 0.0 param['redden'] = 0 obj = Galaxy(param) return obj def defineCCD(iccd, config_file): with open(config_file, "r") as stream: try: config = yaml.safe_load(stream) # for key, value in config.items(): # print (key + " : " + str(value)) except yaml.YAMLError as exc: print(exc) chip = Chip(chipID=iccd, config=config) chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) focal_plane = FocalPlane(chip_list=[iccd]) chip.img.wcs = focal_plane.getTanWCS( 192.8595, 27.1283, -113.4333*galsim.degrees, chip.pix_scale) return chip def defineFilt(chip): filter_param = FilterParam() filter_id, filter_type = chip.getChipFilter() filt = Filter( filter_id=filter_id, filter_type=filter_type, filter_param=filter_param, ccd_bandpass=chip.effCurve) return filt class imagingModule_coverage(unittest.TestCase): def __init__(self, methodName='runTest'): super(imagingModule_coverage, self).__init__(methodName) self.dataPath = os.path.join( os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_msc_sim/csst_fz_msc') self.iccd = 8 def test_imaging(self): config_file = os.path.join(self.dataPath, 'config_test.yaml') chip = defineCCD(self.iccd, config_file) bandpass = defineFilt(chip) filt = defineFilt(chip) print(chip.chipID) print(chip.cen_pix_x, chip.cen_pix_y) catalog = Catalog() obj = catalog._load(os.path.join( self.dataPath, 'galaxies_C6_bundle000287.h5')) sed_data = catalog.load_sed(obj) norm_filt = catalog.load_norm_filt(obj) obj_sed, obj_mag, obj_flux = catalog.convert_sed( mag=obj.param["mag_use_normal"], sed=sed_data, target_filt=filt, norm_filt=norm_filt, ) pupil_area = np.pi * (0.5 * 2.)**2 exptime = 150. # getElectronFluxFilt(filt, tel, exptime) nphotons_tot = obj_flux*pupil_area * exptime full = integrate_sed_bandpass(sed=obj_sed, bandpass=filt.bandpass_full) print(full, nphotons_tot, obj_mag) for i in range(4): sub = integrate_sed_bandpass( sed=obj_sed, bandpass=filt.bandpass_sub_list[i]) ratio = sub / full nphotons = ratio * nphotons_tot disk = galsim.Sersic( n=obj.param['disk_sersic_idx'], half_light_radius=obj.param['hlr_disk'], flux=1.0) disk_shape = galsim.Shear( g1=obj.param['e1_disk'], g2=obj.param['e2_disk']) disk = disk.shear(disk_shape) gal_temp = disk gal_temp = gal_temp.withFlux(nphotons) psf = galsim.Gaussian(sigma=0.1, flux=1) gal_temp = galsim.Convolve(psf, gal_temp) if i == 0: gal = gal_temp else: gal = gal + gal_temp print(gal) self.assertTrue(gal != None) if __name__ == '__main__': unittest.main()