Loading ObservationSim/Config/Header/ImageHeader.py +2 −2 Original line number Diff line number Diff line Loading @@ -376,7 +376,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec h_prim['OBJ_RA'] = ra h_prim['OBJ_DEC'] = dec obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'} obs_type = {'SCI': '01', 'BIAS': '04', 'DARK': '08', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'DARKPF':'09'} OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0') Loading @@ -386,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec h_prim['EXPTIME'] = exptime # Define file types file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'} file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF','DARKPF':'DARK'} h_prim['FILETYPE'] = file_type[im_type] co = coord.SkyCoord(ra, dec, unit='deg') Loading ObservationSim/MockObject/FlatLED.py 0 → 100755 +350 −0 Original line number Diff line number Diff line import galsim import os, sys import numpy as np from astropy.io import fits from scipy.interpolate import griddata import math import astropy.constants as cons from astropy.table import Table from ObservationSim.MockObject.SpecDisperser import SpecDisperser import time from scipy import interpolate from ObservationSim.MockObject.MockObject import MockObject # from ObservationSim.Straylight import calculateSkyMap_split_g # flatDir = '/Volumes/EAGET/LED_FLAT/' LED_name = ['LED1', 'LED2', 'LED3', 'LED4', 'LED5', 'LED6', 'LED7', 'LED8', 'LED9', 'LED10', 'LED11', 'LED12', 'LED13', 'LED14'] cwaves_name = {'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', 'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050', 'LED13': '340', 'LED14': '365'} cwaves = {'LED1': 2750, 'LED2': 3100, 'LED3': 4300, 'LED4': 5050, 'LED5': 5250, 'LED6': 5900, 'LED7': 6700, 'LED8': 7600, 'LED9': 8800, 'LED10': 9400, 'LED11': 10500, 'LED12': 15500, 'LED13': 3400, 'LED14': 3650} cwaves_fwhm = {'LED1': 110, 'LED2': 120, 'LED3': 200, 'LED4': 300, 'LED5': 300, 'LED6': 130, 'LED7': 210, 'LED8': 260, 'LED9': 400, 'LED10': 370, 'LED11': 500, 'LED12': 1400, 'LED13': 90, 'LED14': 100} # LED_QE = {'LED1': 0.3, 'LED2': 0.4, 'LED13': 0.5, 'LED14': 0.5, 'LED10': 0.4} # e-/ms fluxLED = {'LED1': 0.16478729, 'LED2': 0.084220931, 'LED3': 2.263360617, 'LED4': 2.190623489, 'LED5': 0.703504768, 'LED6': 0.446117963, 'LED7': 0.647122098, 'LED8': 0.922313442, 'LED9': 0.987278143, 'LED10': 2.043989167, 'LED11': 0.612571429, 'LED12': 1.228915663, 'LED13': 0.17029384, 'LED14': 0.27842925} mirro_eff = {'GU':0.61, 'GV':0.8, 'GI':0.8} class FlatLED(object): def __init__(self, chip,filt, flatDir = '/Users/zhangxin/Work/SlitlessSim/csst_sls_calibration/flat_field_cube/models/', logger=None): # self.led_type_list = led_type_list self.flatDir = flatDir self.filt = filt self.chip = chip self.logger = logger ### ### return LED flat, e/s ### def getLEDImage(self, led_type='LED1'): # cwave = cwaves[led_type] flat = fits.open(self.flatDir + 'model_' + cwaves_name[led_type] + 'nm.fits') xlen = flat[0].header['NAXIS1'] ylen = 601 x = np.linspace(0, self.chip.npix_x * 6, xlen) y = np.linspace(0, self.chip.npix_y * 5, ylen) xx, yy = np.meshgrid(x, y) a1 = flat[0].data[0:ylen, 0:xlen] # z = np.sin((xx+yy+xx**2+yy**2)) # fInterp = interp2d(xx, yy, z, kind='linear') X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None])) Z_ = a1.flatten() n_x = np.arange(0, self.chip.npix_x * 6, 1) n_y = np.arange(0, self.chip.npix_y * 5, 1) M, N = np.meshgrid(n_x, n_y) i = self.chip.rowID - 1 j = self.chip.colID - 1 U = griddata(X_, Z_, ( M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)], N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]), method='cubic') U = U/np.mean(U) flatImage = U*fluxLED[led_type]*1000 return flatImage def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]): if len(led_type_list) > len(exp_t_list): return np.ones([self.chip.npix_y,self.chip.npix_x]) ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x]) for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] unitFlatImg = self.getLEDImage(led_type=led_type) led_wave = cwaves[led_type] led_fwhm = cwaves_fwhm[led_type] led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm) speci = interpolate.interp1d(led_spec['WAVELENGTH'], led_spec['FLUX']) w_list = np.arange(self.filt.blue_limit, self.filt.red_limit, 0.5) #A f_spec = speci(w_list) ccd_bp = self.chip._getChipEffCurve(self.chip.filter_type) ccd_eff = ccd_bp.__call__(w_list / 10.) filt_bp = self.filt.filter_bandpass fil_eff = filt_bp.__call__(w_list / 10.) t_spec = np.trapz(f_spec*ccd_eff*fil_eff, w_list) print(i, np.mean(unitFlatImg), t_spec, exp_t) unitFlatImg = unitFlatImg * t_spec ledFlat = ledFlat+unitFlatImg*exp_t return ledFlat def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]): if len(led_type_list) != len(exp_t_list): return np.ones([self.chip.npix_y,self.chip.npix_x]) ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x]) for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] unitFlatImg = self.getLEDImage(led_type=led_type) ledFlat_ = unitFlatImg*exp_t ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type] ledFlat_.astype(np.float32) led_wave = cwaves[led_type] led_fwhm = cwaves_fwhm[led_type] led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm) ledspec_map = self.calculateLEDSpec( skyMap=ledFlat_, blueLimit=self.filt.blue_limit, redLimit=self.filt.red_limit, conf=self.chip.sls_conf, pixelSize=self.chip.pix_scale, isAlongY=0, flat_cube=self.chip.flat_cube, led_spec=led_spec) ledFlat = ledFlat + ledspec_map return ledFlat def drawObj_LEDFlat(self, led_type_list=['LED1'], exp_t_list=[0.1]): if self.chip.survey_type == "photometric": return self.drawObj_LEDFlat_img(led_type_list=led_type_list, exp_t_list=exp_t_list) elif self.chip.survey_type == "spectroscopic": return self.drawObj_LEDFlat_slitless(led_type_list=led_type_list, exp_t_list=exp_t_list) def gaussian1d_profile_led(self, xc=5050, fwhm=300): sigma = fwhm/2.355 x_radii = int(5*sigma + 1) xlist = np.arange(xc-x_radii, xc+x_radii, 0.5) xlist_ = np.zeros(len(xlist) + 2) xlist_[1:-1] = xlist xlist_[0] = 2550 xlist_[-1] = 10000 data = np.exp((-(xlist-xc)*(xlist-xc))/(2*sigma*sigma))/(np.sqrt(2*math.pi)*sigma) data_ = np.zeros(len(xlist) + 2) data_[1:-1] = data return Table(np.array([xlist_.astype(np.float32), data_.astype(np.float32)]).T, names=('WAVELENGTH', 'FLUX')) def calculateLEDSpec(self, skyMap=None, blueLimit=4200, redLimit=6500, conf=[''], pixelSize=0.074, isAlongY=0, split_pos=3685, flat_cube=None, led_spec=None): conf1 = conf[0] conf2 = conf[0] if np.size(conf) == 2: conf2 = conf[1] skyImg = galsim.Image(skyMap, xmin=0, ymin=0) tbstart = blueLimit tbend = redLimit fimg = np.zeros_like(skyMap) fImg = galsim.Image(fimg) spec = led_spec if isAlongY == 0: directParm = 0 if isAlongY == 1: directParm = 1 if split_pos >= skyImg.array.shape[directParm]: skyImg1 = galsim.Image(skyImg.array) origin1 = [0, 0] # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, # full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend, # origin=origin1, # conf=conf1) # sdp.compute_spec_orders() y_len = skyMap.shape[0] x_len = skyMap.shape[1] delt_x = 100 delt_y = 100 sub_y_start_arr = np.arange(0, y_len, delt_y) sub_y_end_arr = sub_y_start_arr + delt_y sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len) sub_x_start_arr = np.arange(0, x_len, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf2, flat_cube=flat_cube, ignoreBeam=['D', 'E']) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] else: # sdp.compute_spec_orders() y_len = skyMap.shape[0] x_len = skyMap.shape[1] delt_x = 500 delt_y = y_len sub_y_start_arr = np.arange(0, y_len, delt_y) sub_y_end_arr = sub_y_start_arr + delt_y sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len) delt_x = split_pos - 0 sub_x_start_arr = np.arange(0, split_pos, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf1, flat_cube=flat_cube) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] T2 = time.time() print('time: %s ms' % ((T2 - T1) * 1000)) delt_x = x_len - split_pos sub_x_start_arr = np.arange(split_pos, x_len, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf2, flat_cube=flat_cube) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] T2 = time.time() print('time: %s ms' % ((T2 - T1) * 1000)) if isAlongY == 1: fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0) else: fimg = fImg.array fimg = fimg * pixelSize * pixelSize return fimg ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +1 −1 Original line number Diff line number Diff line Loading @@ -239,7 +239,7 @@ class SpecDisperser(object): # else: # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] status = disperse.disperse_grism_object(self.thumb_img, status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), flat_index[nonz], yfrac_beam[nonz], sensitivity_beam[nonz], modelf, x0, Loading ObservationSim/MockObject/__init__.py +1 −0 Original line number Diff line number Diff line Loading @@ -4,5 +4,6 @@ from .CatalogBase import CatalogBase from .Quasar import Quasar from .Star import Star from .Stamp import Stamp from .FlatLED import FlatLED # from .SkybackgroundMap import * # from .CosmicRay import CosmicRay ObservationSim/ObservationSim.py +115 −7 Original line number Diff line number Diff line Loading @@ -18,6 +18,7 @@ from ObservationSim.Straylight import calculateSkyMap_split_g from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS from ObservationSim._util import get_shear_field, makeSubDir_PointingList from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position from ObservationSim.MockObject import FlatLED class Observation(object): def __init__(self, config, Catalog, work_dir=None, data_dir=None): Loading Loading @@ -202,8 +203,8 @@ class Observation(object): obj = self.cat.objs[j] # (DEBUG) # if obj.getMagFilter(filt)>20: # continue if obj.getMagFilter(filt)>20: continue # load and convert SED; also caculate object's magnitude in all CSST bands try: Loading Loading @@ -412,6 +413,106 @@ class Observation(object): chip_output.Log_info("check running:2: 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) )) def run_one_chip_LEDlight(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, cat_dir=None, sed_dir=None): # # Get WCS for the focal plane # if wcs_fp == None: # wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale) # Create chip Image chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) # chip.img.wcs = wcs_fp led_obj = FlatLED(chip, filt) led_flat = led_obj.drawObj_LEDFlat(led_type_list=self.config["obs_setting"]["LED_TYPE"], exp_t_list=self.config["obs_setting"]["LED_TIME"]) chip.img = chip.img + galsim.ImageF(led_flat) # whether to output zero, dark, flat calibration images. sky_map = np.zeros_like(chip.img.array) self.config["ins_effects"]["shutter_effect"] = False chip.img = chip.addEffects( config=self.config, img=chip.img, chip_output=chip_output, filt=filt, ra_cen=pointing.ra, dec_cen=pointing.dec, img_rot=pointing.img_pa, exptime=self.config["obs_setting"]["exp_time"], pointing_ID=pointing.id, timestamp_obs=pointing.timestamp, pointing_type=pointing.pointing_type, sky_map=sky_map, tel=self.tel, logger=chip_output.logger) if pointing.pointing_type == 'MS': datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") h_prim = generatePrimaryHeader( xlen=chip.npix_x, ylen=chip.npix_y, pointNum=str(pointing.id), ra=pointing.ra, dec=pointing.dec, pixel_scale=chip.pix_scale, date=date_obs, time_obs=time_obs, exptime=self.config["obs_setting"]["exp_time"], im_type='DARKPF', sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], chip_name=str(chip.chipID).rjust(2, '0')) h_ext = generateExtensionHeader( chip=chip, xlen=chip.npix_x, ylen=chip.npix_y, ra=pointing.ra, dec=pointing.dec, pa=pointing.img_pa.deg, gain=chip.gain, readout=chip.read_noise, dark=chip.dark_noise, saturation=90000, pixel_scale=chip.pix_scale, pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, extName='SCI', timestamp=pointing.timestamp, exptime=self.config["obs_setting"]["exp_time"], readoutTime=chip.readout_time) chip.img = galsim.Image(chip.img.array, dtype=np.uint16) hdu1 = fits.PrimaryHDU(header=h_prim) hdu1.add_checksum() hdu1.header.comments['CHECKSUM'] = 'HDU checksum' hdu1.header.comments['DATASUM'] = 'data unit checksum' hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) hdu2.add_checksum() hdu2.header.comments['XTENSION'] = 'extension type' hdu2.header.comments['CHECKSUM'] = 'HDU checksum' hdu2.header.comments['DATASUM'] = 'data unit checksum' hdu1 = fits.HDUList([hdu1, hdu2]) fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits') hdu1.writeto(fname, output_verify='ignore', overwrite=True) # chip_output.Log_info("# objects that are too bright %d out of %d" % (bright_obj, self.nobj)) # chip_output.Log_info("# objects that are too dim %d out of %d" % (dim_obj, self.nobj)) # chip_output.Log_info("# objects that are missed %d out of %d" % (missed_obj, self.nobj)) del chip.img chip_output.Log_info("check running:2: 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))) def runExposure_MPI_PointingList(self, pointing_list,chips=None, use_mpi=False): if use_mpi: comm = MPI.COMM_WORLD Loading @@ -434,6 +535,7 @@ class Observation(object): run_chips.append(chip) run_filts.append(filt) for ipoint in range(len(pointing_list)): for ichip in range(nchips_per_fp): i = ipoint*nchips_per_fp + ichip Loading Loading @@ -461,6 +563,12 @@ class Observation(object): subdir=sub_img_dir, prefix=prefix) chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid)) if self.config["obs_setting"]["survey_type"] == "LED": self.run_one_chip_LEDlight(chip=chip, filt=filt, chip_output=chip_output, pointing=pointing) else: self.run_one_chip( chip=chip, filt=filt, Loading Loading
ObservationSim/Config/Header/ImageHeader.py +2 −2 Original line number Diff line number Diff line Loading @@ -376,7 +376,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec h_prim['OBJ_RA'] = ra h_prim['OBJ_DEC'] = dec obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'} obs_type = {'SCI': '01', 'BIAS': '04', 'DARK': '08', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'DARKPF':'09'} OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0') Loading @@ -386,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec h_prim['EXPTIME'] = exptime # Define file types file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'} file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF','DARKPF':'DARK'} h_prim['FILETYPE'] = file_type[im_type] co = coord.SkyCoord(ra, dec, unit='deg') Loading
ObservationSim/MockObject/FlatLED.py 0 → 100755 +350 −0 Original line number Diff line number Diff line import galsim import os, sys import numpy as np from astropy.io import fits from scipy.interpolate import griddata import math import astropy.constants as cons from astropy.table import Table from ObservationSim.MockObject.SpecDisperser import SpecDisperser import time from scipy import interpolate from ObservationSim.MockObject.MockObject import MockObject # from ObservationSim.Straylight import calculateSkyMap_split_g # flatDir = '/Volumes/EAGET/LED_FLAT/' LED_name = ['LED1', 'LED2', 'LED3', 'LED4', 'LED5', 'LED6', 'LED7', 'LED8', 'LED9', 'LED10', 'LED11', 'LED12', 'LED13', 'LED14'] cwaves_name = {'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', 'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050', 'LED13': '340', 'LED14': '365'} cwaves = {'LED1': 2750, 'LED2': 3100, 'LED3': 4300, 'LED4': 5050, 'LED5': 5250, 'LED6': 5900, 'LED7': 6700, 'LED8': 7600, 'LED9': 8800, 'LED10': 9400, 'LED11': 10500, 'LED12': 15500, 'LED13': 3400, 'LED14': 3650} cwaves_fwhm = {'LED1': 110, 'LED2': 120, 'LED3': 200, 'LED4': 300, 'LED5': 300, 'LED6': 130, 'LED7': 210, 'LED8': 260, 'LED9': 400, 'LED10': 370, 'LED11': 500, 'LED12': 1400, 'LED13': 90, 'LED14': 100} # LED_QE = {'LED1': 0.3, 'LED2': 0.4, 'LED13': 0.5, 'LED14': 0.5, 'LED10': 0.4} # e-/ms fluxLED = {'LED1': 0.16478729, 'LED2': 0.084220931, 'LED3': 2.263360617, 'LED4': 2.190623489, 'LED5': 0.703504768, 'LED6': 0.446117963, 'LED7': 0.647122098, 'LED8': 0.922313442, 'LED9': 0.987278143, 'LED10': 2.043989167, 'LED11': 0.612571429, 'LED12': 1.228915663, 'LED13': 0.17029384, 'LED14': 0.27842925} mirro_eff = {'GU':0.61, 'GV':0.8, 'GI':0.8} class FlatLED(object): def __init__(self, chip,filt, flatDir = '/Users/zhangxin/Work/SlitlessSim/csst_sls_calibration/flat_field_cube/models/', logger=None): # self.led_type_list = led_type_list self.flatDir = flatDir self.filt = filt self.chip = chip self.logger = logger ### ### return LED flat, e/s ### def getLEDImage(self, led_type='LED1'): # cwave = cwaves[led_type] flat = fits.open(self.flatDir + 'model_' + cwaves_name[led_type] + 'nm.fits') xlen = flat[0].header['NAXIS1'] ylen = 601 x = np.linspace(0, self.chip.npix_x * 6, xlen) y = np.linspace(0, self.chip.npix_y * 5, ylen) xx, yy = np.meshgrid(x, y) a1 = flat[0].data[0:ylen, 0:xlen] # z = np.sin((xx+yy+xx**2+yy**2)) # fInterp = interp2d(xx, yy, z, kind='linear') X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None])) Z_ = a1.flatten() n_x = np.arange(0, self.chip.npix_x * 6, 1) n_y = np.arange(0, self.chip.npix_y * 5, 1) M, N = np.meshgrid(n_x, n_y) i = self.chip.rowID - 1 j = self.chip.colID - 1 U = griddata(X_, Z_, ( M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)], N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]), method='cubic') U = U/np.mean(U) flatImage = U*fluxLED[led_type]*1000 return flatImage def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]): if len(led_type_list) > len(exp_t_list): return np.ones([self.chip.npix_y,self.chip.npix_x]) ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x]) for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] unitFlatImg = self.getLEDImage(led_type=led_type) led_wave = cwaves[led_type] led_fwhm = cwaves_fwhm[led_type] led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm) speci = interpolate.interp1d(led_spec['WAVELENGTH'], led_spec['FLUX']) w_list = np.arange(self.filt.blue_limit, self.filt.red_limit, 0.5) #A f_spec = speci(w_list) ccd_bp = self.chip._getChipEffCurve(self.chip.filter_type) ccd_eff = ccd_bp.__call__(w_list / 10.) filt_bp = self.filt.filter_bandpass fil_eff = filt_bp.__call__(w_list / 10.) t_spec = np.trapz(f_spec*ccd_eff*fil_eff, w_list) print(i, np.mean(unitFlatImg), t_spec, exp_t) unitFlatImg = unitFlatImg * t_spec ledFlat = ledFlat+unitFlatImg*exp_t return ledFlat def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]): if len(led_type_list) != len(exp_t_list): return np.ones([self.chip.npix_y,self.chip.npix_x]) ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x]) for i in np.arange(len(led_type_list)): led_type = led_type_list[i] exp_t = exp_t_list[i] unitFlatImg = self.getLEDImage(led_type=led_type) ledFlat_ = unitFlatImg*exp_t ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type] ledFlat_.astype(np.float32) led_wave = cwaves[led_type] led_fwhm = cwaves_fwhm[led_type] led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm) ledspec_map = self.calculateLEDSpec( skyMap=ledFlat_, blueLimit=self.filt.blue_limit, redLimit=self.filt.red_limit, conf=self.chip.sls_conf, pixelSize=self.chip.pix_scale, isAlongY=0, flat_cube=self.chip.flat_cube, led_spec=led_spec) ledFlat = ledFlat + ledspec_map return ledFlat def drawObj_LEDFlat(self, led_type_list=['LED1'], exp_t_list=[0.1]): if self.chip.survey_type == "photometric": return self.drawObj_LEDFlat_img(led_type_list=led_type_list, exp_t_list=exp_t_list) elif self.chip.survey_type == "spectroscopic": return self.drawObj_LEDFlat_slitless(led_type_list=led_type_list, exp_t_list=exp_t_list) def gaussian1d_profile_led(self, xc=5050, fwhm=300): sigma = fwhm/2.355 x_radii = int(5*sigma + 1) xlist = np.arange(xc-x_radii, xc+x_radii, 0.5) xlist_ = np.zeros(len(xlist) + 2) xlist_[1:-1] = xlist xlist_[0] = 2550 xlist_[-1] = 10000 data = np.exp((-(xlist-xc)*(xlist-xc))/(2*sigma*sigma))/(np.sqrt(2*math.pi)*sigma) data_ = np.zeros(len(xlist) + 2) data_[1:-1] = data return Table(np.array([xlist_.astype(np.float32), data_.astype(np.float32)]).T, names=('WAVELENGTH', 'FLUX')) def calculateLEDSpec(self, skyMap=None, blueLimit=4200, redLimit=6500, conf=[''], pixelSize=0.074, isAlongY=0, split_pos=3685, flat_cube=None, led_spec=None): conf1 = conf[0] conf2 = conf[0] if np.size(conf) == 2: conf2 = conf[1] skyImg = galsim.Image(skyMap, xmin=0, ymin=0) tbstart = blueLimit tbend = redLimit fimg = np.zeros_like(skyMap) fImg = galsim.Image(fimg) spec = led_spec if isAlongY == 0: directParm = 0 if isAlongY == 1: directParm = 1 if split_pos >= skyImg.array.shape[directParm]: skyImg1 = galsim.Image(skyImg.array) origin1 = [0, 0] # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, # full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend, # origin=origin1, # conf=conf1) # sdp.compute_spec_orders() y_len = skyMap.shape[0] x_len = skyMap.shape[1] delt_x = 100 delt_y = 100 sub_y_start_arr = np.arange(0, y_len, delt_y) sub_y_end_arr = sub_y_start_arr + delt_y sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len) sub_x_start_arr = np.arange(0, x_len, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf2, flat_cube=flat_cube, ignoreBeam=['D', 'E']) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] else: # sdp.compute_spec_orders() y_len = skyMap.shape[0] x_len = skyMap.shape[1] delt_x = 500 delt_y = y_len sub_y_start_arr = np.arange(0, y_len, delt_y) sub_y_end_arr = sub_y_start_arr + delt_y sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len) delt_x = split_pos - 0 sub_x_start_arr = np.arange(0, split_pos, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf1, flat_cube=flat_cube) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] T2 = time.time() print('time: %s ms' % ((T2 - T1) * 1000)) delt_x = x_len - split_pos sub_x_start_arr = np.arange(split_pos, x_len, delt_x) sub_x_end_arr = sub_x_start_arr + delt_x sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len) for i, k1 in enumerate(sub_y_start_arr): sub_y_s = k1 sub_y_e = sub_y_end_arr[i] sub_y_center = (sub_y_s + sub_y_e) / 2. for j, k2 in enumerate(sub_x_start_arr): sub_x_s = k2 sub_x_e = sub_x_end_arr[j] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e]) origin_sub = [sub_y_s, sub_x_s] sub_x_center = (sub_x_s + sub_x_e) / 2. sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub, tar_spec=spec, band_start=tbstart, band_end=tbend, conf=conf2, flat_cube=flat_cube) spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] origin_order_x = v[1] origin_order_y = v[2] ssImg = galsim.ImageF(img_s) ssImg.setOrigin(origin_order_x, origin_order_y) bounds = ssImg.bounds & fImg.bounds if bounds.area() == 0: continue fImg[bounds] = fImg[bounds] + ssImg[bounds] T2 = time.time() print('time: %s ms' % ((T2 - T1) * 1000)) if isAlongY == 1: fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0) else: fimg = fImg.array fimg = fimg * pixelSize * pixelSize return fimg
ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +1 −1 Original line number Diff line number Diff line Loading @@ -239,7 +239,7 @@ class SpecDisperser(object): # else: # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] status = disperse.disperse_grism_object(self.thumb_img, status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), flat_index[nonz], yfrac_beam[nonz], sensitivity_beam[nonz], modelf, x0, Loading
ObservationSim/MockObject/__init__.py +1 −0 Original line number Diff line number Diff line Loading @@ -4,5 +4,6 @@ from .CatalogBase import CatalogBase from .Quasar import Quasar from .Star import Star from .Stamp import Stamp from .FlatLED import FlatLED # from .SkybackgroundMap import * # from .CosmicRay import CosmicRay
ObservationSim/ObservationSim.py +115 −7 Original line number Diff line number Diff line Loading @@ -18,6 +18,7 @@ from ObservationSim.Straylight import calculateSkyMap_split_g from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS from ObservationSim._util import get_shear_field, makeSubDir_PointingList from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position from ObservationSim.MockObject import FlatLED class Observation(object): def __init__(self, config, Catalog, work_dir=None, data_dir=None): Loading Loading @@ -202,8 +203,8 @@ class Observation(object): obj = self.cat.objs[j] # (DEBUG) # if obj.getMagFilter(filt)>20: # continue if obj.getMagFilter(filt)>20: continue # load and convert SED; also caculate object's magnitude in all CSST bands try: Loading Loading @@ -412,6 +413,106 @@ class Observation(object): chip_output.Log_info("check running:2: 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) )) def run_one_chip_LEDlight(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, cat_dir=None, sed_dir=None): # # Get WCS for the focal plane # if wcs_fp == None: # wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale) # Create chip Image chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) # chip.img.wcs = wcs_fp led_obj = FlatLED(chip, filt) led_flat = led_obj.drawObj_LEDFlat(led_type_list=self.config["obs_setting"]["LED_TYPE"], exp_t_list=self.config["obs_setting"]["LED_TIME"]) chip.img = chip.img + galsim.ImageF(led_flat) # whether to output zero, dark, flat calibration images. sky_map = np.zeros_like(chip.img.array) self.config["ins_effects"]["shutter_effect"] = False chip.img = chip.addEffects( config=self.config, img=chip.img, chip_output=chip_output, filt=filt, ra_cen=pointing.ra, dec_cen=pointing.dec, img_rot=pointing.img_pa, exptime=self.config["obs_setting"]["exp_time"], pointing_ID=pointing.id, timestamp_obs=pointing.timestamp, pointing_type=pointing.pointing_type, sky_map=sky_map, tel=self.tel, logger=chip_output.logger) if pointing.pointing_type == 'MS': datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") h_prim = generatePrimaryHeader( xlen=chip.npix_x, ylen=chip.npix_y, pointNum=str(pointing.id), ra=pointing.ra, dec=pointing.dec, pixel_scale=chip.pix_scale, date=date_obs, time_obs=time_obs, exptime=self.config["obs_setting"]["exp_time"], im_type='DARKPF', sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], chip_name=str(chip.chipID).rjust(2, '0')) h_ext = generateExtensionHeader( chip=chip, xlen=chip.npix_x, ylen=chip.npix_y, ra=pointing.ra, dec=pointing.dec, pa=pointing.img_pa.deg, gain=chip.gain, readout=chip.read_noise, dark=chip.dark_noise, saturation=90000, pixel_scale=chip.pix_scale, pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, extName='SCI', timestamp=pointing.timestamp, exptime=self.config["obs_setting"]["exp_time"], readoutTime=chip.readout_time) chip.img = galsim.Image(chip.img.array, dtype=np.uint16) hdu1 = fits.PrimaryHDU(header=h_prim) hdu1.add_checksum() hdu1.header.comments['CHECKSUM'] = 'HDU checksum' hdu1.header.comments['DATASUM'] = 'data unit checksum' hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) hdu2.add_checksum() hdu2.header.comments['XTENSION'] = 'extension type' hdu2.header.comments['CHECKSUM'] = 'HDU checksum' hdu2.header.comments['DATASUM'] = 'data unit checksum' hdu1 = fits.HDUList([hdu1, hdu2]) fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits') hdu1.writeto(fname, output_verify='ignore', overwrite=True) # chip_output.Log_info("# objects that are too bright %d out of %d" % (bright_obj, self.nobj)) # chip_output.Log_info("# objects that are too dim %d out of %d" % (dim_obj, self.nobj)) # chip_output.Log_info("# objects that are missed %d out of %d" % (missed_obj, self.nobj)) del chip.img chip_output.Log_info("check running:2: 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))) def runExposure_MPI_PointingList(self, pointing_list,chips=None, use_mpi=False): if use_mpi: comm = MPI.COMM_WORLD Loading @@ -434,6 +535,7 @@ class Observation(object): run_chips.append(chip) run_filts.append(filt) for ipoint in range(len(pointing_list)): for ichip in range(nchips_per_fp): i = ipoint*nchips_per_fp + ichip Loading Loading @@ -461,6 +563,12 @@ class Observation(object): subdir=sub_img_dir, prefix=prefix) chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid)) if self.config["obs_setting"]["survey_type"] == "LED": self.run_one_chip_LEDlight(chip=chip, filt=filt, chip_output=chip_output, pointing=pointing) else: self.run_one_chip( chip=chip, filt=filt, Loading