From 88aea0d6162878a9917188e1213b64fece6a3450 Mon Sep 17 00:00:00 2001
From: fangyuedong <fangyuedong@qxslab.cn>
Date: Tue, 2 Apr 2024 10:42:48 +0800
Subject: [PATCH] split sim steps into multiple files
---
ObservationSim/ObservationSim.py | 2 +-
ObservationSim/SimSteps.py | 547 ------------------
ObservationSim/sim_steps/__init__.py | 35 ++
ObservationSim/sim_steps/add_LED_flat.py | 19 +
.../sim_steps/add_brighter_fatter_CTE.py | 42 ++
ObservationSim/sim_steps/add_cosmic_rays.py | 11 +
ObservationSim/sim_steps/add_objects.py | 196 +++++++
ObservationSim/sim_steps/add_pattern_noise.py | 72 +++
.../sim_steps/add_sky_background.py | 101 ++++
ObservationSim/sim_steps/prepare_headers.py | 38 ++
ObservationSim/sim_steps/readout_output.py | 61 ++
11 files changed, 576 insertions(+), 548 deletions(-)
delete mode 100644 ObservationSim/SimSteps.py
create mode 100644 ObservationSim/sim_steps/__init__.py
create mode 100644 ObservationSim/sim_steps/add_LED_flat.py
create mode 100644 ObservationSim/sim_steps/add_brighter_fatter_CTE.py
create mode 100644 ObservationSim/sim_steps/add_cosmic_rays.py
create mode 100644 ObservationSim/sim_steps/add_objects.py
create mode 100644 ObservationSim/sim_steps/add_pattern_noise.py
create mode 100644 ObservationSim/sim_steps/add_sky_background.py
create mode 100644 ObservationSim/sim_steps/prepare_headers.py
create mode 100644 ObservationSim/sim_steps/readout_output.py
diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py
index f790224..278b560 100755
--- a/ObservationSim/ObservationSim.py
+++ b/ObservationSim/ObservationSim.py
@@ -14,7 +14,7 @@ from ObservationSim.Instrument.Chip import Effects
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim._util import makeSubDir_PointingList
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
-from ObservationSim.SimSteps import SimSteps, SIM_STEP_TYPES
+from ObservationSim.sim_steps import SimSteps, SIM_STEP_TYPES
class Observation(object):
def __init__(self, config, Catalog, work_dir=None, data_dir=None):
diff --git a/ObservationSim/SimSteps.py b/ObservationSim/SimSteps.py
deleted file mode 100644
index 5300c52..0000000
--- a/ObservationSim/SimSteps.py
+++ /dev/null
@@ -1,547 +0,0 @@
-import os
-import galsim
-import traceback
-import gc
-import psutil
-import numpy as np
-from astropy.io import fits
-from datetime import datetime
-from numpy.random import Generator, PCG64
-
-from ObservationSim._util import get_shear_field
-from ObservationSim.Straylight import calculateSkyMap_split_g
-from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
-from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS
-from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
-from ObservationSim.Instrument.Chip import Effects
-from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim
-
-from ObservationSim.MockObject import FlatLED
-from ObservationSim.Instrument.FilterParam import FilterParam
-
-class SimSteps:
- def __init__(self, overall_config, chip_output, all_filters):
- self.overall_config = overall_config
- self.chip_output = chip_output
- self.all_filters = all_filters
-
- def prepare_headers(self, chip, pointing):
-
- self.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,
- time_pt = pointing.timestamp,
- exptime=pointing.exp_time,
- im_type=pointing.pointing_type,
- sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
- sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
- project_cycle=self.overall_config["project_cycle"],
- run_counter=self.overall_config["run_counter"],
- chip_name=str(chip.chipID).rjust(2, '0'))
- self.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=pointing.pointing_type,
- timestamp = pointing.timestamp,
- exptime = pointing.exp_time,
- readoutTime = chip.readout_time)
- return self.h_prim, self.h_ext
-
- def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param):
-
- if not hasattr(self, 'h_ext'):
- _, _ = self.prepare_headers(chip=chip, pointing=pointing)
- chip_wcs = galsim.FitsWCS(header = self.h_ext)
-
- expTime = obs_param["exptime"]
- skyback_level = obs_param["flat_level"]
-
- filter_param = FilterParam()
- sky_level_filt = obs_param["flat_level_filt"]
- norm_scaler = skyback_level/expTime / filter_param.param[sky_level_filt][5]
-
- flat_normal = np.ones_like(chip.img.array)
- if obs_param["flat_fielding"] == True:
- flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
- if obs_param["shutter_effect"] == True:
- flat_normal = flat_normal * chip.shutter_img
- flat_normal = np.array(flat_normal, dtype='float32')
-
-
- if chip.survey_type == "photometric":
- sky_map = flat_normal * np.ones_like(chip.img.array) * norm_scaler * filter_param.param[chip.filter_type][5] / tel.pupil_area * expTime
- elif chip.survey_type == "spectroscopic":
- # flat_normal = np.ones_like(chip.img.array)
- if obs_param["flat_fielding"] == True:
-
- flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
- if obs_param["shutter_effect"] == True:
-
- flat_normal = flat_normal * chip.shutter_img
- 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 * expTime
-
- chip.img += sky_map
- return chip, filt, tel, pointing
-
- def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
- if not hasattr(self, 'h_ext'):
- _, _ = self.prepare_headers(chip=chip, pointing=pointing)
- chip_wcs = galsim.FitsWCS(header = self.h_ext)
-
- if "flat_level" not in obs_param or "flat_level_filt" not in obs_param:
- chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
- else:
- if obs_param.get('flat_level') is None or obs_param.get('flat_level_filt')is None:
- chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
- else:
- chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)
-
- # chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
-
- return chip, filt, tel, pointing
-
-
-
- def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
- flat_normal = np.ones_like(chip.img.array)
- if obs_param["flat_fielding"] == True:
- flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
- if obs_param["shutter_effect"] == True:
- flat_normal = flat_normal * chip.shutter_img
- flat_normal = np.array(flat_normal, dtype='float32')
-
- if obs_param["enable_straylight_model"]:
- # Filter.sky_background, Filter.zodical_spec will be updated
- filt.setFilterStrayLightPixel(
- jtime = pointing.jdt,
- sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]),
- pointing_radec = np.array([pointing.ra,pointing.dec]),
- sun_pos = np.array([pointing.sun_x, pointing.sun_y, pointing.sun_z]))
- self.chip_output.Log_info("================================================")
- self.chip_output.Log_info("sky background + stray light pixel flux value: %.5f"%(filt.sky_background))
-
- if chip.survey_type == "photometric":
- sky_map = filt.getSkyNoise(exptime = obs_param["exptime"])
- sky_map = sky_map * np.ones_like(chip.img.array) * flat_normal
- sky_map = galsim.Image(array=sky_map)
- else:
- # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits')
- 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,
- zoldial_spec = filt.zodical_spec)
- sky_map = (sky_map + filt.sky_background)*obs_param["exptime"]
-
- # sky_map = sky_map * tel.pupil_area * obs_param["exptime"]
- chip.img += sky_map
- return chip, filt, tel, pointing
-
- def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
-
- if not hasattr(self, 'h_ext'):
- _, _ = self.prepare_headers(chip=chip, pointing=pointing)
- chip_wcs = galsim.FitsWCS(header = self.h_ext)
- pf_map = np.zeros_like(chip.img.array)
- if obs_param["LED_TYPE"] is not None:
- if len(obs_param["LED_TYPE"]) != 0:
- print("LED OPEN--------")
-
- led_obj = FlatLED(chip, filt)
-
- led_flat = led_obj.drawObj_LEDFlat(led_type_list=obs_param["LED_TYPE"], exp_t_list=obs_param["LED_TIME"])
- pf_map = led_flat
- chip.img = chip.img + led_flat
- return chip, filt, tel, pointing
-
- def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
-
- # Prepare output file(s) for this chip
- self.chip_output.create_output_file()
-
- # Prepare the PSF model
- if self.overall_config["psf_setting"]["psf_model"] == "Gauss":
- psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"])
- elif self.overall_config["psf_setting"]["psf_model"] == "Interp":
- if chip.survey_type == "spectroscopic":
- psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"])
- else:
- psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"])
- else:
- self.chip_output.Log_error("unrecognized PSF model type!!", flush=True)
-
- # Apply field distortion model
- if obs_param["field_dist"] == True:
- fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg)
- else:
- fd_model = None
-
- # Update limiting magnitudes for all filters based on the exposure time
- # Get the filter which will be used for magnitude cut
- for ifilt in range(len(self.all_filters)):
- temp_filter = self.all_filters[ifilt]
- 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():
- cut_filter = temp_filter
-
- # Read in shear values from configuration file if the constant shear type is used
- if self.overall_config["shear_setting"]["shear_type"] == "constant":
- 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))
-
- # Get chip WCS
- if not hasattr(self, 'h_ext'):
- _, _ = self.prepare_headers(chip=chip, pointing=pointing)
- chip_wcs = galsim.FitsWCS(header = self.h_ext)
-
- # Loop over objects
- nobj = len(catalog.objs)
- missed_obj = 0
- bright_obj = 0
- dim_obj = 0
- for j in range(nobj):
- # # [DEBUG] [TODO]
- # if j >= 10:
- # break
- obj = catalog.objs[j]
-
- # load and convert SED; also caculate object's magnitude in all CSST bands
- try:
- sed_data = catalog.load_sed(obj)
- norm_filt = catalog.load_norm_filt(obj)
- obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = catalog.convert_sed(
- mag=obj.param["mag_use_normal"],
- sed=sed_data,
- target_filt=filt,
- norm_filt=norm_filt,
- )
- _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = catalog.convert_sed(
- mag=obj.param["mag_use_normal"],
- sed=sed_data,
- target_filt=cut_filter,
- norm_filt=norm_filt,
- )
- except Exception as e:
- traceback.print_exc()
- self.chip_output.Log_error(e)
- continue
-
- # [TODO] Testing
- # self.chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()]))
-
- # Exclude very bright/dim objects (for now)
- if cut_filter.is_too_bright(
- mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()],
- 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()]))
- bright_obj += 1
- obj.unload_SED()
- continue
- if filt.is_too_dim(
- mag=obj.getMagFilter(filt),
- 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)))
- dim_obj += 1
- obj.unload_SED()
- continue
-
- # Get corresponding shear values
- if self.overall_config["shear_setting"]["shear_type"] == "constant":
- if obj.type == 'star':
- obj.g1, obj.g2 = 0., 0.
- else:
- # Figure out shear fields from overall configuration shear setting
- obj.g1, obj.g2 = g1_field, g2_field
- elif self.overall_config["shear_setting"]["shear_type"] == "catalog":
- pass
- else:
- self.chip_output.Log_error("Unknown shear input")
- raise ValueError("Unknown shear input")
-
- # 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)
-
- # [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
- # 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:
- 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_error("Objected missed: %s"%(obj.id))
- missed_obj += 1
- obj.unload_SED()
- continue
-
- # Draw object & update output catalog
- try:
- if self.overall_config["run_option"]["out_cat_only"]:
- 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)
- pos_shear = 0.
- elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]:
- isUpdated, pos_shear = obj.drawObj_multiband(
- tel=tel,
- pos_img=pos_img,
- psf_model=psf_model,
- bandpass_list=filt.bandpass_sub_list,
- filt=filt,
- chip=chip,
- g1=obj.g1,
- g2=obj.g2,
- exptime=pointing.exp_time,
- fd_shear=fd_shear)
-
- elif chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"]:
- isUpdated, pos_shear = obj.drawObj_slitless(
- tel=tel,
- pos_img=pos_img,
- psf_model=psf_model,
- bandpass_list=filt.bandpass_sub_list,
- filt=filt,
- chip=chip,
- g1=obj.g1,
- g2=obj.g2,
- exptime=pointing.exp_time,
- normFilter=norm_filt,
- fd_shear=fd_shear)
-
- if isUpdated == 1:
- # TODO: add up stats
- self.chip_output.cat_add_obj(obj, pos_img, pos_shear)
- pass
- elif isUpdated == 0:
- missed_obj += 1
- self.chip_output.Log_error("Objected missed: %s"%(obj.id))
- else:
- self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id))
- continue
- except Exception as e:
- traceback.print_exc()
- self.chip_output.Log_error(e)
-
- # Unload SED:
- obj.unload_SED()
- del obj
- gc.collect()
- del psf_model
- 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("# objects that are too bright %d out of %d"%(bright_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)
- flat_normal = np.ones_like(chip.img.array)
- if obs_param["flat_fielding"] == True:
- flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
- if obs_param["shutter_effect"] == True:
- flat_normal = flat_normal * chip.shutter_img
- flat_normal = np.array(flat_normal, dtype='float32')
- chip.img *= flat_normal
- del flat_normal
-
- return chip, filt, tel, pointing
-
- def add_cosmic_rays(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Adding Cosmic-Ray")
- chip.img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(
- img=chip.img,
- chip=chip,
- exptime=pointing.exp_time,
- seed=self.overall_config["random_seeds"]["seed_CR"]+pointing.id*30+chip.chipID)
- # [TODO] output cosmic ray image
- return chip, filt, tel, pointing
-
- def apply_PRNU(self, chip, filt, tel, pointing, catalog, obs_param):
- chip.img *= chip.prnu_img
- return chip, filt, tel, pointing
-
- def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
- # Add dark current & Poisson noise
- InputDark = False
- if obs_param["add_dark"] == True:
- if InputDark:
- chip.img = chip_utils.add_inputdark(img=chip.img, chip=chip, exptime=pointing.exp_time)
- else:
- chip.img, _ = chip_utils.add_poisson(img=chip.img, chip=chip, exptime=pointing.exp_time, poisson_noise=chip.poisson_noise)
- else:
- chip.img, _ = chip_utils.add_poisson(img=chip.img, chip=self, exptime=pointing.exp_time, poisson_noise=chip.poisson_noise, dark_noise=0.)
- return chip, filt, tel, pointing
-
- def add_brighter_fatter(self, chip, filt, tel, pointing, catalog, obs_param):
- chip.img = chip_utils.add_brighter_fatter(img=chip.img)
- return chip, filt, tel, pointing
-
- def add_detector_defects(self, chip, filt, tel, pointing, catalog, obs_param):
- # Add Hot Pixels or/and Dead Pixels
- rgbadpix = Generator(PCG64(int(self.overall_config["random_seeds"]["seed_defective"]+chip.chipID)))
- badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
- chip.img = Effects.DefectivePixels(
- chip.img,
- IfHotPix=obs_param["hot_pixels"],
- IfDeadPix=obs_param["dead_pixels"],
- fraction=badfraction,
- seed=self.overall_config["random_seeds"]["seed_defective"]+chip.chipID, biaslevel=0)
- # Apply Bad columns
- if obs_param["bad_columns"] == True:
- chip.img = Effects.BadColumns(chip.img, seed=self.overall_config["random_seeds"]["seed_badcolumns"], chipid=chip.chipID)
- return chip, filt, tel, pointing
-
- def add_nonlinearity(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Applying Non-Linearity on the chip image")
- chip.img = Effects.NonLinearity(GSImage=chip.img, beta1=5.e-7, beta2=0)
- return chip, filt, tel, pointing
-
- def add_blooming(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Applying CCD Saturation & Blooming")
- chip.img = Effects.SaturBloom(GSImage=chip.img, nsect_x=1, nsect_y=1, fullwell=int(chip.full_well))
- return chip, filt, tel, pointing
-
- def apply_CTE(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Apply CTE Effect")
- ### 2*8 -> 1*16 img-layout
- img = chip_utils.formatOutput(GSImage=chip.img)
- chip.nsecy = 1
- chip.nsecx = 16
-
- img_arr = img.array
- ny, nx = img_arr.shape
- dx = int(nx/chip.nsecx)
- dy = int(ny/chip.nsecy)
- newimg = galsim.Image(nx, int(ny+chip.overscan_y), init_value=0)
- for ichannel in range(16):
- print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing.id, chip.chipID, ichannel+1))
- noverscan, nsp, nmax = chip.overscan_y, 3, 10
- beta, w, c = 0.478, 84700, 0
- t = np.array([0.74, 7.7, 37],dtype=np.float32)
- rho_trap = np.array([0.6, 1.6, 1.4],dtype=np.float32)
- trap_seeds = np.array([0, 1000, 10000],dtype=np.int32) + ichannel + chip.chipID*16
- release_seed = 50 + ichannel + pointing.id*30 + chip.chipID*16
- newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
- newimg.wcs = img.wcs
- del img
- img = newimg
-
- ### 1*16 -> 2*8 img-layout
- chip.img = chip_utils.formatRevert(GSImage=img)
- chip.nsecy = 2
- chip.nsecx = 8
-
- # [TODO] make overscan_y == 0
- chip.overscan_y = 0
- return chip, filt, tel, pointing
-
- def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info("Apply pre/over-scan")
- chip.img = chip_utils.AddPreScan(GSImage=chip.img,
- pre1=chip.prescan_x,
- pre2=chip.prescan_y,
- over1=chip.overscan_x,
- over2=chip.overscan_y)
- return chip, filt, tel, pointing
-
- def add_bias(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Adding Bias level and 16-channel non-uniformity")
- if obs_param["bias_16channel"] == True:
- chip.img = Effects.AddBiasNonUniform16(chip.img,
- bias_level=float(chip.bias_level),
- nsecy = chip.nsecy, nsecx=chip.nsecx,
- seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID)
- elif obs_param["bias_16channel"] == False:
- chip.img += self.bias_level
- return chip, filt, tel, pointing
-
- def add_readout_noise(self, chip, filt, tel, pointing, catalog, obs_param):
- seed = int(self.overall_config["random_seeds"]["seed_readout"]) + pointing.id*30 + chip.chipID
- rng_readout = galsim.BaseDeviate(seed)
- readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=chip.read_noise)
- chip.img.addNoise(readout_noise)
- return chip, filt, tel, pointing
-
- def apply_gain(self, chip, filt, tel, pointing, catalog, obs_param):
- self.chip_output.Log_info(" Applying Gain")
- if obs_param["gain_16channel"] == True:
- chip.img, chip.gain_channel = Effects.ApplyGainNonUniform16(
- chip.img, gain=chip.gain,
- nsecy = chip.nsecy, nsecx=chip.nsecx,
- seed=self.overall_config["random_seeds"]["seed_gainNonUniform"]+chip.chipID)
- elif obs_param["gain_16channel"] == False:
- chip.img /= chip.gain
- return chip, filt, tel, pointing
-
- def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param):
-
- if obs_param["format_output"] == True:
- self.chip_output.Log_info(" Apply 1*16 format")
- chip.img = chip_utils.formatOutput(GSImage=chip.img)
- chip.nsecy = 1
- chip.nsecx = 16
-
- chip.img.array[chip.img.array > 65535] = 65535
- chip.img.replaceNegative(replace_value=0)
- chip.img.quantize()
- chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
- hdu1 = fits.PrimaryHDU(header=self.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=self.h_ext)
- hdu2.add_checksum()
- hdu2.header.comments['XTENSION'] = 'extension type'
- hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
- hdu2.header.comments['DATASUM'] = 'data unit checksum'
- hdu2.header.comments["XTENSION"] = "image extension"
- hdu1 = fits.HDUList([hdu1, hdu2])
- fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits')
- hdu1.writeto(fname, output_verify='ignore', overwrite=True)
- return chip, filt, tel, pointing
-
-SIM_STEP_TYPES = {
- "scie_obs": "add_objects",
- "sky_background": "add_sky_background",
- "cosmic_rays": "add_cosmic_rays",
- "PRNU_effect": "apply_PRNU",
- "poisson_and_dark": "add_poisson_and_dark",
- "bright_fatter": "add_brighter_fatter",
- "detector_defects": "add_detector_defects",
- "nonlinearity": "add_nonlinearity",
- "blooming": "add_blooming",
- "CTE_effect": "apply_CTE",
- "prescan_overscan": "add_prescan_overscan",
- "bias": "add_bias",
- "readout_noise": "add_readout_noise",
- "gain": "apply_gain",
- "quantization_and_output": "quantization_and_output",
- "led_calib_model":"add_LED_Flat",
- # "sky_flatField":"add_sky_flat_calibration",
-}
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/__init__.py b/ObservationSim/sim_steps/__init__.py
new file mode 100644
index 0000000..e5ccadf
--- /dev/null
+++ b/ObservationSim/sim_steps/__init__.py
@@ -0,0 +1,35 @@
+import os
+
+class SimSteps:
+ def __init__(self, overall_config, chip_output, all_filters):
+ self.overall_config = overall_config
+ self.chip_output = chip_output
+ self.all_filters = all_filters
+
+ from .prepare_headers import prepare_headers
+ from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background
+ from .add_objects import add_objects
+ from .add_cosmic_rays import add_cosmic_rays
+ from .add_pattern_noise import apply_PRNU, add_poisson_and_dark, add_detector_defects, add_nonlinearity, add_blooming, add_bias
+ from .add_brighter_fatter_CTE import add_brighter_fatter, apply_CTE
+ from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output
+
+SIM_STEP_TYPES = {
+ "scie_obs": "add_objects",
+ "sky_background": "add_sky_background",
+ "cosmic_rays": "add_cosmic_rays",
+ "PRNU_effect": "apply_PRNU",
+ "poisson_and_dark": "add_poisson_and_dark",
+ "bright_fatter": "add_brighter_fatter",
+ "detector_defects": "add_detector_defects",
+ "nonlinearity": "add_nonlinearity",
+ "blooming": "add_blooming",
+ "CTE_effect": "apply_CTE",
+ "prescan_overscan": "add_prescan_overscan",
+ "bias": "add_bias",
+ "readout_noise": "add_readout_noise",
+ "gain": "apply_gain",
+ "quantization_and_output": "quantization_and_output",
+ "led_calib_model":"add_LED_Flat",
+ "sky_flatField":"add_sky_flat_calibration",
+}
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_LED_flat.py b/ObservationSim/sim_steps/add_LED_flat.py
new file mode 100644
index 0000000..b8a5ac2
--- /dev/null
+++ b/ObservationSim/sim_steps/add_LED_flat.py
@@ -0,0 +1,19 @@
+import numpy as np
+from ObservationSim.MockObject import FlatLED
+
+def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
+
+ if not hasattr(self, 'h_ext'):
+ _, _ = self.prepare_headers(chip=chip, pointing=pointing)
+ chip_wcs = galsim.FitsWCS(header = self.h_ext)
+ pf_map = np.zeros_like(chip.img.array)
+ if obs_param["LED_TYPE"] is not None:
+ if len(obs_param["LED_TYPE"]) != 0:
+ print("LED OPEN--------")
+
+ led_obj = FlatLED(chip, filt)
+
+ led_flat = led_obj.drawObj_LEDFlat(led_type_list=obs_param["LED_TYPE"], exp_t_list=obs_param["LED_TIME"])
+ pf_map = led_flat
+ chip.img = chip.img + led_flat
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_brighter_fatter_CTE.py b/ObservationSim/sim_steps/add_brighter_fatter_CTE.py
new file mode 100644
index 0000000..ec8a31a
--- /dev/null
+++ b/ObservationSim/sim_steps/add_brighter_fatter_CTE.py
@@ -0,0 +1,42 @@
+import numpy as np
+import galsim
+from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
+from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim
+
+def add_brighter_fatter(self, chip, filt, tel, pointing, catalog, obs_param):
+ chip.img = chip_utils.add_brighter_fatter(img=chip.img)
+ return chip, filt, tel, pointing
+
+def apply_CTE(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Apply CTE Effect")
+ ### 2*8 -> 1*16 img-layout
+ img = chip_utils.formatOutput(GSImage=chip.img)
+ chip.nsecy = 1
+ chip.nsecx = 16
+
+ img_arr = img.array
+ ny, nx = img_arr.shape
+ dx = int(nx/chip.nsecx)
+ dy = int(ny/chip.nsecy)
+ newimg = galsim.Image(nx, int(ny+chip.overscan_y), init_value=0)
+ for ichannel in range(16):
+ print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing.id, chip.chipID, ichannel+1))
+ noverscan, nsp, nmax = chip.overscan_y, 3, 10
+ beta, w, c = 0.478, 84700, 0
+ t = np.array([0.74, 7.7, 37],dtype=np.float32)
+ rho_trap = np.array([0.6, 1.6, 1.4],dtype=np.float32)
+ trap_seeds = np.array([0, 1000, 10000],dtype=np.int32) + ichannel + chip.chipID*16
+ release_seed = 50 + ichannel + pointing.id*30 + chip.chipID*16
+ newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
+ newimg.wcs = img.wcs
+ del img
+ img = newimg
+
+ ### 1*16 -> 2*8 img-layout
+ chip.img = chip_utils.formatRevert(GSImage=img)
+ chip.nsecy = 2
+ chip.nsecx = 8
+
+ # [TODO] make overscan_y == 0
+ chip.overscan_y = 0
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_cosmic_rays.py b/ObservationSim/sim_steps/add_cosmic_rays.py
new file mode 100644
index 0000000..48f8bd8
--- /dev/null
+++ b/ObservationSim/sim_steps/add_cosmic_rays.py
@@ -0,0 +1,11 @@
+from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
+
+def add_cosmic_rays(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Adding Cosmic-Ray")
+ chip.img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(
+ img=chip.img,
+ chip=chip,
+ exptime=pointing.exp_time,
+ seed=self.overall_config["random_seeds"]["seed_CR"]+pointing.id*30+chip.chipID)
+ # [TODO] output cosmic ray image
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_objects.py b/ObservationSim/sim_steps/add_objects.py
new file mode 100644
index 0000000..05d96cf
--- /dev/null
+++ b/ObservationSim/sim_steps/add_objects.py
@@ -0,0 +1,196 @@
+import os
+import gc
+import psutil
+import numpy as np
+import galsim
+from ObservationSim._util import get_shear_field
+from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS
+
+def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
+
+ # Prepare output file(s) for this chip
+ self.chip_output.create_output_file()
+
+ # Prepare the PSF model
+ if self.overall_config["psf_setting"]["psf_model"] == "Gauss":
+ psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"])
+ elif self.overall_config["psf_setting"]["psf_model"] == "Interp":
+ if chip.survey_type == "spectroscopic":
+ psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"])
+ else:
+ psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"])
+ else:
+ self.chip_output.Log_error("unrecognized PSF model type!!", flush=True)
+
+ # Apply field distortion model
+ if obs_param["field_dist"] == True:
+ fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg)
+ else:
+ fd_model = None
+
+ # Update limiting magnitudes for all filters based on the exposure time
+ # Get the filter which will be used for magnitude cut
+ for ifilt in range(len(self.all_filters)):
+ temp_filter = self.all_filters[ifilt]
+ 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():
+ cut_filter = temp_filter
+
+ # Read in shear values from configuration file if the constant shear type is used
+ if self.overall_config["shear_setting"]["shear_type"] == "constant":
+ 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))
+
+ # Get chip WCS
+ if not hasattr(self, 'h_ext'):
+ _, _ = self.prepare_headers(chip=chip, pointing=pointing)
+ chip_wcs = galsim.FitsWCS(header = self.h_ext)
+
+ # Loop over objects
+ nobj = len(catalog.objs)
+ missed_obj = 0
+ bright_obj = 0
+ dim_obj = 0
+ for j in range(nobj):
+ # # [DEBUG] [TODO]
+ # if j >= 10:
+ # break
+ obj = catalog.objs[j]
+
+ # load and convert SED; also caculate object's magnitude in all CSST bands
+ try:
+ sed_data = catalog.load_sed(obj)
+ norm_filt = catalog.load_norm_filt(obj)
+ obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = catalog.convert_sed(
+ mag=obj.param["mag_use_normal"],
+ sed=sed_data,
+ target_filt=filt,
+ norm_filt=norm_filt,
+ )
+ _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = catalog.convert_sed(
+ mag=obj.param["mag_use_normal"],
+ sed=sed_data,
+ target_filt=cut_filter,
+ norm_filt=norm_filt,
+ )
+ except Exception as e:
+ traceback.print_exc()
+ self.chip_output.Log_error(e)
+ continue
+
+ # [TODO] Testing
+ # self.chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()]))
+
+ # Exclude very bright/dim objects (for now)
+ if cut_filter.is_too_bright(
+ mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()],
+ 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()]))
+ bright_obj += 1
+ obj.unload_SED()
+ continue
+ if filt.is_too_dim(
+ mag=obj.getMagFilter(filt),
+ 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)))
+ dim_obj += 1
+ obj.unload_SED()
+ continue
+
+ # Get corresponding shear values
+ if self.overall_config["shear_setting"]["shear_type"] == "constant":
+ if obj.type == 'star':
+ obj.g1, obj.g2 = 0., 0.
+ else:
+ # Figure out shear fields from overall configuration shear setting
+ obj.g1, obj.g2 = g1_field, g2_field
+ elif self.overall_config["shear_setting"]["shear_type"] == "catalog":
+ pass
+ else:
+ self.chip_output.Log_error("Unknown shear input")
+ raise ValueError("Unknown shear input")
+
+ # 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)
+
+ # [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
+ # 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:
+ 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_error("Objected missed: %s"%(obj.id))
+ missed_obj += 1
+ obj.unload_SED()
+ continue
+
+ # Draw object & update output catalog
+ try:
+ if self.overall_config["run_option"]["out_cat_only"]:
+ 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)
+ pos_shear = 0.
+ elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]:
+ isUpdated, pos_shear = obj.drawObj_multiband(
+ tel=tel,
+ pos_img=pos_img,
+ psf_model=psf_model,
+ bandpass_list=filt.bandpass_sub_list,
+ filt=filt,
+ chip=chip,
+ g1=obj.g1,
+ g2=obj.g2,
+ exptime=pointing.exp_time,
+ fd_shear=fd_shear)
+
+ elif chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"]:
+ isUpdated, pos_shear = obj.drawObj_slitless(
+ tel=tel,
+ pos_img=pos_img,
+ psf_model=psf_model,
+ bandpass_list=filt.bandpass_sub_list,
+ filt=filt,
+ chip=chip,
+ g1=obj.g1,
+ g2=obj.g2,
+ exptime=pointing.exp_time,
+ normFilter=norm_filt,
+ fd_shear=fd_shear)
+
+ if isUpdated == 1:
+ # TODO: add up stats
+ self.chip_output.cat_add_obj(obj, pos_img, pos_shear)
+ pass
+ elif isUpdated == 0:
+ missed_obj += 1
+ self.chip_output.Log_error("Objected missed: %s"%(obj.id))
+ else:
+ self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id))
+ continue
+ except Exception as e:
+ traceback.print_exc()
+ self.chip_output.Log_error(e)
+
+ # Unload SED:
+ obj.unload_SED()
+ del obj
+ gc.collect()
+ del psf_model
+ 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("# objects that are too bright %d out of %d"%(bright_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)
+ flat_normal = np.ones_like(chip.img.array)
+ if obs_param["flat_fielding"] == True:
+ flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
+ if obs_param["shutter_effect"] == True:
+ flat_normal = flat_normal * chip.shutter_img
+ flat_normal = np.array(flat_normal, dtype='float32')
+ chip.img *= flat_normal
+ del flat_normal
+
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_pattern_noise.py b/ObservationSim/sim_steps/add_pattern_noise.py
new file mode 100644
index 0000000..ceb07d4
--- /dev/null
+++ b/ObservationSim/sim_steps/add_pattern_noise.py
@@ -0,0 +1,72 @@
+from numpy.random import Generator, PCG64
+from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
+from ObservationSim.Instrument.Chip import Effects
+
+def apply_PRNU(self, chip, filt, tel, pointing, catalog, obs_param):
+ chip.img *= chip.prnu_img
+ return chip, filt, tel, pointing
+
+def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
+ # Add dark current & Poisson noise
+ InputDark = False
+ if obs_param["add_dark"] == True:
+ if InputDark:
+ chip.img = chip_utils.add_inputdark(img=chip.img,
+ chip=chip,
+ exptime=pointing.exp_time)
+ else:
+ chip.img, _ = chip_utils.add_poisson(img=chip.img,
+ chip=chip,
+ exptime=pointing.exp_time,
+ poisson_noise=chip.poisson_noise)
+ else:
+ chip.img, _ = chip_utils.add_poisson(img=chip.img,
+ chip=self,
+ exptime=pointing.exp_time,
+ poisson_noise=chip.poisson_noise,
+ dark_noise=0.)
+ return chip, filt, tel, pointing
+
+def add_detector_defects(self, chip, filt, tel, pointing, catalog, obs_param):
+ # Add Hot Pixels or/and Dead Pixels
+ rgbadpix = Generator(PCG64(int(self.overall_config["random_seeds"]["seed_defective"]+chip.chipID)))
+ badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
+ chip.img = Effects.DefectivePixels(
+ chip.img,
+ IfHotPix=obs_param["hot_pixels"],
+ IfDeadPix=obs_param["dead_pixels"],
+ fraction=badfraction,
+ seed=self.overall_config["random_seeds"]["seed_defective"]+chip.chipID, biaslevel=0)
+ # Apply Bad columns
+ if obs_param["bad_columns"] == True:
+ chip.img = Effects.BadColumns(chip.img,
+ seed=self.overall_config["random_seeds"]["seed_badcolumns"],
+ chipid=chip.chipID)
+ return chip, filt, tel, pointing
+
+def add_nonlinearity(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Applying Non-Linearity on the chip image")
+ chip.img = Effects.NonLinearity(GSImage=chip.img,
+ beta1=5.e-7,
+ beta2=0)
+ return chip, filt, tel, pointing
+
+def add_blooming(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Applying CCD Saturation & Blooming")
+ chip.img = Effects.SaturBloom(GSImage=chip.img,
+ nsect_x=1,
+ nsect_y=1,
+ fullwell=int(chip.full_well))
+ return chip, filt, tel, pointing
+
+def add_bias(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Adding Bias level and 16-channel non-uniformity")
+ if obs_param["bias_16channel"] == True:
+ chip.img = Effects.AddBiasNonUniform16(chip.img,
+ bias_level=float(chip.bias_level),
+ nsecy = chip.nsecy,
+ nsecx=chip.nsecx,
+ seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID)
+ elif obs_param["bias_16channel"] == False:
+ chip.img += self.bias_level
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/add_sky_background.py b/ObservationSim/sim_steps/add_sky_background.py
new file mode 100644
index 0000000..fc50878
--- /dev/null
+++ b/ObservationSim/sim_steps/add_sky_background.py
@@ -0,0 +1,101 @@
+import numpy as np
+import galsim
+from ObservationSim.Straylight import calculateSkyMap_split_g
+
+def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
+ flat_normal = np.ones_like(chip.img.array)
+ if obs_param["flat_fielding"] == True:
+ flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
+ if obs_param["shutter_effect"] == True:
+ flat_normal = flat_normal * chip.shutter_img
+ flat_normal = np.array(flat_normal, dtype='float32')
+
+ if obs_param["enable_straylight_model"]:
+ # Filter.sky_background, Filter.zodical_spec will be updated
+ filt.setFilterStrayLightPixel(
+ jtime = pointing.jdt,
+ sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]),
+ pointing_radec = np.array([pointing.ra,pointing.dec]),
+ sun_pos = np.array([pointing.sun_x, pointing.sun_y, pointing.sun_z]))
+ self.chip_output.Log_info("================================================")
+ self.chip_output.Log_info("sky background + stray light pixel flux value: %.5f"%(filt.sky_background))
+
+ if chip.survey_type == "photometric":
+ sky_map = filt.getSkyNoise(exptime = obs_param["exptime"])
+ sky_map = sky_map * np.ones_like(chip.img.array) * flat_normal
+ sky_map = galsim.Image(array=sky_map)
+ else:
+ # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits')
+ 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,
+ zoldial_spec = filt.zodical_spec)
+ sky_map = (sky_map + filt.sky_background)*obs_param["exptime"]
+
+ # sky_map = sky_map * tel.pupil_area * obs_param["exptime"]
+ chip.img += sky_map
+ return chip, filt, tel, pointing
+
+def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param):
+
+ if not hasattr(self, 'h_ext'):
+ _, _ = self.prepare_headers(chip=chip, pointing=pointing)
+ chip_wcs = galsim.FitsWCS(header = self.h_ext)
+
+ expTime = obs_param["exptime"]
+ skyback_level = obs_param["flat_level"]
+
+ filter_param = FilterParam()
+ sky_level_filt = obs_param["flat_level_filt"]
+ norm_scaler = skyback_level/expTime / filter_param.param[sky_level_filt][5]
+
+ flat_normal = np.ones_like(chip.img.array)
+ if obs_param["flat_fielding"] == True:
+ flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
+ if obs_param["shutter_effect"] == True:
+ flat_normal = flat_normal * chip.shutter_img
+ flat_normal = np.array(flat_normal, dtype='float32')
+
+
+ if chip.survey_type == "photometric":
+ sky_map = flat_normal * np.ones_like(chip.img.array) * norm_scaler * filter_param.param[chip.filter_type][5] / tel.pupil_area * expTime
+ elif chip.survey_type == "spectroscopic":
+ # flat_normal = np.ones_like(chip.img.array)
+ if obs_param["flat_fielding"] == True:
+
+ flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
+ if obs_param["shutter_effect"] == True:
+
+ flat_normal = flat_normal * chip.shutter_img
+ 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 * expTime
+
+ chip.img += sky_map
+ return chip, filt, tel, pointing
+
+def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
+ if not hasattr(self, 'h_ext'):
+ _, _ = self.prepare_headers(chip=chip, pointing=pointing)
+ chip_wcs = galsim.FitsWCS(header = self.h_ext)
+
+ if "flat_level" not in obs_param or "flat_level_filt" not in obs_param:
+ chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
+ else:
+ if obs_param.get('flat_level') is None or obs_param.get('flat_level_filt')is None:
+ chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
+ else:
+ chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)
+ return chip, filt, tel, pointing
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/prepare_headers.py b/ObservationSim/sim_steps/prepare_headers.py
new file mode 100644
index 0000000..935add5
--- /dev/null
+++ b/ObservationSim/sim_steps/prepare_headers.py
@@ -0,0 +1,38 @@
+from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
+
+def prepare_headers(self, chip, pointing):
+ self.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,
+ time_pt = pointing.timestamp,
+ exptime=pointing.exp_time,
+ im_type=pointing.pointing_type,
+ sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
+ sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
+ project_cycle=self.overall_config["project_cycle"],
+ run_counter=self.overall_config["run_counter"],
+ chip_name=str(chip.chipID).rjust(2, '0'))
+ self.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=pointing.pointing_type,
+ timestamp = pointing.timestamp,
+ exptime = pointing.exp_time,
+ readoutTime = chip.readout_time)
+ return self.h_prim, self.h_ext
\ No newline at end of file
diff --git a/ObservationSim/sim_steps/readout_output.py b/ObservationSim/sim_steps/readout_output.py
new file mode 100644
index 0000000..380d445
--- /dev/null
+++ b/ObservationSim/sim_steps/readout_output.py
@@ -0,0 +1,61 @@
+import os
+import galsim
+import numpy as np
+from astropy.io import fits
+from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
+from ObservationSim.Instrument.Chip import Effects
+
+def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info("Apply pre/over-scan")
+ chip.img = chip_utils.AddPreScan(GSImage=chip.img,
+ pre1=chip.prescan_x,
+ pre2=chip.prescan_y,
+ over1=chip.overscan_x,
+ over2=chip.overscan_y)
+ return chip, filt, tel, pointing
+
+def add_readout_noise(self, chip, filt, tel, pointing, catalog, obs_param):
+ seed = int(self.overall_config["random_seeds"]["seed_readout"]) + pointing.id*30 + chip.chipID
+ rng_readout = galsim.BaseDeviate(seed)
+ readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=chip.read_noise)
+ chip.img.addNoise(readout_noise)
+ return chip, filt, tel, pointing
+
+def apply_gain(self, chip, filt, tel, pointing, catalog, obs_param):
+ self.chip_output.Log_info(" Applying Gain")
+ if obs_param["gain_16channel"] == True:
+ chip.img, chip.gain_channel = Effects.ApplyGainNonUniform16(chip.img,
+ gain=chip.gain,
+ nsecy = chip.nsecy,
+ nsecx=chip.nsecx,
+ seed=self.overall_config["random_seeds"]["seed_gainNonUniform"]+chip.chipID)
+ elif obs_param["gain_16channel"] == False:
+ chip.img /= chip.gain
+ return chip, filt, tel, pointing
+
+def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param):
+
+ if obs_param["format_output"] == True:
+ self.chip_output.Log_info(" Apply 1*16 format")
+ chip.img = chip_utils.formatOutput(GSImage=chip.img)
+ chip.nsecy = 1
+ chip.nsecx = 16
+
+ chip.img.array[chip.img.array > 65535] = 65535
+ chip.img.replaceNegative(replace_value=0)
+ chip.img.quantize()
+ chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
+ hdu1 = fits.PrimaryHDU(header=self.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=self.h_ext)
+ hdu2.add_checksum()
+ hdu2.header.comments['XTENSION'] = 'extension type'
+ hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
+ hdu2.header.comments['DATASUM'] = 'data unit checksum'
+ hdu2.header.comments["XTENSION"] = "image extension"
+ hdu1 = fits.HDUList([hdu1, hdu2])
+ fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits')
+ hdu1.writeto(fname, output_verify='ignore', overwrite=True)
+ return chip, filt, tel, pointing
--
GitLab