diff --git a/ObservationSim/Config/Pointing.py b/ObservationSim/Config/Pointing.py index 6c17ea4a414d5ec494fc81e43c667bf47e70f987..dccf91861b85bd48e9e741577d9e4e094ac606cf 100644 --- a/ObservationSim/Config/Pointing.py +++ b/ObservationSim/Config/Pointing.py @@ -59,7 +59,7 @@ class Pointing(object): self.sat_z = float(columns[8]) self.sun_x = float(columns[9]) self.sun_y = float(columns[10]) - self.sun_z = float(columns[1]) + self.sun_z = float(columns[11]) self.sat_vx = float(columns[15]) self.sat_vy = float(columns[16]) self.sat_vz = float(columns[17]) diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 36e3fd6e64bb48037a817138debd697f8fb8c973..2e1c9cab2dfc49323d631acc779bd7e5f391e222 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -270,7 +270,7 @@ class Chip(FocalPlane): noise = self.dark_noise * exptime + self.read_noise**2 return noise - def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, tel=None, logger=None): + def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, post_flash_map=None, tel=None, logger=None): # Set random seeds SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) @@ -309,6 +309,9 @@ class Chip(FocalPlane): if config["output_setting"]["flat_output"] == False: del flat_img + if post_flash_map is not None: + img = img + post_flash_map + # Apply Shutter-effect for one chip if config["ins_effects"]["shutter_effect"] == True: chip_utils.log_info(msg=" Apply shutter effect", logger=self.logger) @@ -320,7 +323,6 @@ class Chip(FocalPlane): shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID)) del shutt_gsimg del shuttimg - # # Add Poisson noise to the resulting images # # (NOTE): this can only applied to the slitless image # # since it dose not use photon shooting to draw stamps diff --git a/ObservationSim/MockObject/FlatLED.py b/ObservationSim/MockObject/FlatLED.py new file mode 100755 index 0000000000000000000000000000000000000000..e21edff0b88532fe87dd7ff59e4d7b84faf10195 --- /dev/null +++ b/ObservationSim/MockObject/FlatLED.py @@ -0,0 +1,350 @@ + + +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 + + + + diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py index 7c8b51efa9f39b07e485f00506a53eaff163f118..d4126f0c457dff809a06912c5c313c3dfda7adf8 100644 --- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +++ b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py @@ -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, diff --git a/ObservationSim/MockObject/__init__.py b/ObservationSim/MockObject/__init__.py index 8faeb8a8c8736b6603c8dc24f8212d82c6741f9e..261107321d91c55eb6424d4dcdbceddfc4957ea3 100755 --- a/ObservationSim/MockObject/__init__.py +++ b/ObservationSim/MockObject/__init__.py @@ -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 diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 9500e473bfa57738c2902be20e335acd0ce7e33f..8d14414a720a4a82be10ca1e062e3b3ce5f80300 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -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): @@ -196,8 +197,8 @@ class Observation(object): for j in range(self.nobj): # (DEBUG) - # if j >= 10: - # break + if j >= 10: + break obj = self.cat.objs[j] @@ -414,6 +415,139 @@ 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_calibration(self, chip, filt, pointing, chip_output, skyback_level = 20000, sky_level_filt = 'g', 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 + pf_map = np.zeros_like(chip.img.array) + if self.config["obs_setting"]["LED_TYPE"] is not None: + if len(self.config["obs_setting"]["LED_TYPE"]) != 0: + print("LED OPEN--------") + + 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"]) + pf_map = led_flat + # whether to output zero, dark, flat calibration images. + expTime = self.config["obs_setting"]["exp_time"] + skybg_unit = self.filter_param.param[sky_level_filt][5] + norm_scaler = skyback_level/expTime/skybg_unit + + if skyback_level == 0: + self.config["ins_effects"]["shutter_effect"] = False + + if chip.survey_type == "photometric": + sky_map = np.ones_like(chip.img.array) * skybg_unit * norm_scaler / self.tel.pupil_area + elif chip.survey_type == "spectroscopic": + flat_normal = np.ones_like(chip.img.array) + if self.config["ins_effects"]["flat_fielding"] == True: + chip_output.Log_info("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding" % chip.chipID) + msg = str(chip.img.bounds) + chip_output.Log_info(msg) + flat_img = Effects.MakeFlatSmooth( + chip.img.bounds, + int(self.config["random_seeds"]["seed_flat"])) + flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array) + if self.config["ins_effects"]["shutter_effect"] == True: + chip_output.Log_info("SLS flat preprocess,CHIP %d : Apply shutter effect" % chip.chipID) + shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735, + dt=1E-3) # shutter effect normalized image for this chip + flat_normal = flat_normal * shuttimg + flat_normal = np.array(flat_normal, dtype='float32') + sky_map = calculateSkyMap_split_g( + skyMap=flat_normal, + blueLimit=filt.blue_limit, + redLimit=filt.red_limit, + conf=chip.sls_conf, + pixelSize=chip.pix_scale, + isAlongY=0, + flat_cube=chip.flat_cube) + sky_map = sky_map * norm_scaler + + 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, + post_flash_map=pf_map, + logger=chip_output.logger) + + 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 @@ -436,6 +570,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 @@ -463,11 +598,19 @@ 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)) - self.run_one_chip( - chip=chip, - filt=filt, - chip_output=chip_output, - pointing=pointing) + if self.config["obs_setting"]["survey_type"] == "CALIBRATION": + self.run_one_chip_calibration(chip=chip, + filt=filt, + chip_output=chip_output, + pointing=pointing, + skyback_level = self.config["obs_setting"]["FLAT_LEVEL"], + sky_level_filt = self.config["obs_setting"]["FLAT_LEVEL_FIL"]) + else: + self.run_one_chip( + chip=chip, + filt=filt, + chip_output=chip_output, + pointing=pointing) chip_output.Log_info("finished running chip#%d..."%(chip.chipID)) for handler in chip_output.logger.handlers[:]: chip_output.logger.removeHandler(handler) diff --git a/config/config_C6_dev.yaml b/config/config_C6_dev.yaml index c1b5ba1a74b2485a6776c07eac44805d80168b83..c7b4442df54c53422d6f618dc6fd25384dc3f198 100644 --- a/config/config_C6_dev.yaml +++ b/config/config_C6_dev.yaml @@ -9,14 +9,10 @@ # Base diretories and naming setup # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent -##<<<<<<< HEAD -##work_dir: "/share/home/zhangxin/CSST_SIM/CSST_new_sim/csst-simulation/" -##======= -work_dir: "/share/home/weichengliang/CSST_git/test_new_sim/outputs/" -##>>>>>>> new_sim -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "testRun2" -project_cycle: 8 +work_dir: "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/" +data_dir: "/Volumes/EAGET/C6_data/inputData/" +run_name: "C6_new_sim_2sq_run1" +project_cycle: 6 run_counter: 1 # Whether to use MPI @@ -70,7 +66,16 @@ obs_setting: # "Spectroscopic": simulate slitless spectroscopic chips only # "FGS": simulate FGS chips only (31-42) # "All": simulate full focal plane - survey_type: "Photometric" + # "CALIBRATION": falt, bias, dark with or without postflash + survey_type: "CALIBRATION" + #"LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null + #'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' + LED_TYPE: ['LED5'] + LED_TIME: [1.] + # unit e- ,flat level + FLAT_LEVEL: 20000 + FLAT_LEVEL_FIL: 'g' # Exposure time [seconds] exp_time: 150. @@ -90,7 +95,7 @@ obs_setting: # if you just want to run default pointing: # - pointing_dir: null # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" + pointing_dir: "/Volumes/EAGET/C6_data/inputData/" pointing_file: "pointing_radec_246.5_40.dat" # Number of calibration pointings @@ -106,10 +111,10 @@ obs_setting: # - give a list of indexes of chips: [ip_1, ip_2...] # - run all chips: null # Note: for all pointings - run_chips: [8] + run_chips: [5] # Whether to enable astrometric modeling - enable_astrometric_model: True + enable_astrometric_model: False # Whether to enable straylight model enable_straylight_model: True @@ -132,7 +137,7 @@ psf_setting: # Which PSF model to use: # "Gauss": simple gaussian profile # "Interp": Interpolated PSF from sampled ray-tracing data - psf_model: "Interp" + psf_model: "Gauss" # PSF size [arcseconds] # radius of 80% energy encircled @@ -164,11 +169,12 @@ ins_effects: # switches # Note: bias_16channel, gain_16channel, and shutter_effect # is currently not applicable to "FGS" observations - field_dist: YES # Whether to add field distortions + field_dist: NO # Whether to add field distortions add_back: YES # Whether to add sky background add_dark: YES # Whether to add dark noise add_readout: YES # Whether to add read-out (Gaussian) noise add_bias: YES # Whether to add bias-level to images + add_prescan: OFF bias_16channel: YES # Whether to add different biases for 16 channels gain_16channel: YES # Whether to make different gains for 16 channels shutter_effect: YES # Whether to add shutter effect