diff --git a/config/obs_config_SCI.yaml b/config/obs_config_SCI.yaml index b80039cfb5fb80f969e8f9c17c582d8fc85deaec..2fc05f96efebd5592f9284eeafeaa91f1adfa3c7 100644 --- a/config/obs_config_SCI.yaml +++ b/config/obs_config_SCI.yaml @@ -9,13 +9,14 @@ ############################################### # Observation type -obs_type: null +obs_type: "WIDE" obs_type_code: "101" obs_id: "00000001" # this setting will only be used if pointing list file is not given # Define list of chips -run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips -# run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips +# run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips +#run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips +run_chips: [17] # Define observation sequence call_sequence: @@ -27,6 +28,11 @@ call_sequence: shutter_effect: YES flat_fielding: YES field_dist: YES + + ghost_SCI: + mag_threshold: 16 + field_dist: YES + # Accumulate fluxes from sky background sky_background: # [Optional]: exposure time of the pointing will be used as default. @@ -83,5 +89,5 @@ call_sequence: gain_16channel: YES # Output the final image quantization_and_output: - format_output: YES + format_output: NO ... diff --git a/observation_sim/mock_objects/Ghost.py b/observation_sim/mock_objects/Ghost.py new file mode 100644 index 0000000000000000000000000000000000000000..14392bc635f53fa181f901555cba2602b933f4e8 --- /dev/null +++ b/observation_sim/mock_objects/Ghost.py @@ -0,0 +1,199 @@ +""" +Ghost image calculation (single-class refactor) + +This module keeps the original behavior but wraps everything into ONE class: `Ghost`. + +All lengths are in meters. +""" + +from __future__ import annotations + +from typing import Dict, Tuple + +import numpy as np +import galsim + + +class Ghost: + """ + Single-class ghost image model. + + Usage + ----- + ghost = Ghost() # uses built-in defaults + result = ghost.calculate(detector_id=13, x0=0.0, y0=0.0) + + The returned dict contains: + - "ghosts": ndarray shape (3, 4) with rows [x, y, radius, relative_energy] + - "angle_x_deg": float + - "angle_y_deg": float + - "wave_name": str (band name) + """ + + # ----------- Default optical constants (meters) ----------- + DEFAULT_ZP = 4.593 + DEFAULT_RP = 0.16367 + DEFAULT_D1 = 0.005 + DEFAULT_D2 = 0.008 + DEFAULT_PIXEL = 1.0e-5 # kept for completeness; unused by formula + + def __init__( + self, + *, + Zp: float = DEFAULT_ZP, + Rp: float = DEFAULT_RP, + D1: float = DEFAULT_D1, + D2: float = DEFAULT_D2, + pixel: float = DEFAULT_PIXEL, + ) -> None: + self.Zp = float(Zp) + self.Rp = float(Rp) + self.D1 = float(D1) + self.D2 = float(D2) + self.pixel = float(pixel) + + # ---- Wave table (from original script) ---- + # Stored as: wave_num -> (wave_name, wave_length, nd, tf1, rf1, tf2, rf2, TD, RD) + self._waves: Dict[int, Tuple[str, float, float, float, float, float, float, float, float]] = { + 0: ("nuv", 2.9e-7, 1.495, 0.825, 0.0093, 0.825, 0.0093, 0.555, 0.15), + 1: ("u", 3.6e-7, 1.473, 0.954, 0.0138, 0.954, 0.0138, 0.56, 0.15), + 2: ("g", 4.8e-7, 1.463, 0.98, 0.0197, 0.98, 0.0197, 0.8, 0.2), + 3: ("r", 6.2e-7, 1.458, 0.98, 0.015, 0.98, 0.015, 0.91, 0.09), + 4: ("i", 7.5e-7, 1.454, 0.985, 0.01, 0.985, 0.01, 0.889, 0.11), + 5: ("z", 9.1e-7, 1.452, 0.99, 0.0085, 0.99, 0.0085, 0.6, 0.1), + 6: ("y", 9.7e-7, 1.451, 0.99, 0.0035, 0.99, 0.0035, 0.36, 0.1), + } + + # ---- Detector -> wave_num mapping (from original script) ---- + self._detector_to_wave_num: Dict[int, int] = { + 6: 6, + 7: 4, + 8: 2, + 9: 3, + 11: 5, + 12: 0, + 13: 0, + 14: 1, + 15: 6, + 16: 6, + 17: 1, + 18: 0, + 19: 0, + 20: 5, + 22: 3, + 23: 2, + 24: 4, + 25: 6, + } + + # ----------------------------- + # Helpers + # ----------------------------- + + def _get_wave_tuple(self, detector_id: int) -> Tuple[str, float, float, float, float, float, float, float, float]: + """ + Return wave properties tuple for a detector_id. + Raises KeyError with a friendly message if unknown. + """ + try: + wave_num = self._detector_to_wave_num[detector_id] + except KeyError as e: + known = ", ".join(map(str, sorted(self._detector_to_wave_num.keys()))) + raise KeyError(f"Unknown detector_id={detector_id}. Known ids: {known}") from e + + return self._waves[wave_num] + + # ----------------------------- + # Public API + # ----------------------------- + + def calculate(self, detector_id: int, x0: float, y0: float) -> Dict[str, object]: + """ + Calculate ghost images for a given detector and focal-plane coordinate. + + Parameters + ---------- + detector_id : int + Detector identifier. + x0, y0 : float + Absolute focal-plane coordinates of the image point [m]. + + Returns + ------- + dict with keys: + ghosts: ndarray (3, 4) rows [ghost_x, ghost_y, ghost_radius, relative_energy] + angle_x_deg: float + angle_y_deg: float + wave_name: str + """ + wave_name, _wave_length, nd, tf1, rf1, tf2, rf2, TD, RD = self._get_wave_tuple(detector_id) + + # Slope defined by pupil geometry + u = self.Rp / self.Zp + + # Ghost radii (physical lengths) + r1 = 2.0 * self.D1 * u / nd + r2 = 2.0 * self.D2 * u + r3 = 2.0 * (self.D1 / nd + self.D2) * u + + # Relative energy terms + m1 = rf1 * rf2 + m2 = rf2 * RD + m3 = rf1 * RD * (tf2 ** 2) + + # NOTE: Original script uses a fixed +0.313 x-offset; preserve behavior. + denom = self.Zp - self.D1 * (1.0 - 1.0 / nd) + angle_x = (float(x0) + 0.313) / denom + angle_y = float(y0) / denom + + # Angles inside the filter + af_x = angle_x / nd + af_y = angle_y / nd + + # Lateral shifts on the focal plane + xd1 = np.tan(af_x) * 2.0 * self.D1 + yd1 = np.tan(af_y) * 2.0 * self.D1 + xd2 = np.tan(angle_x) * 2.0 * self.D2 + yd2 = np.tan(angle_y) * 2.0 * self.D2 + + # Combined shifts + xd3 = xd1 + xd2 + yd3 = yd1 + yd2 + + # Ghost positions + x1, y1 = float(x0) + xd1, float(y0) + yd1 + x2, y2 = float(x0) + xd2, float(y0) + yd2 + x3, y3 = float(x0) + xd3, float(y0) + yd3 + + ghosts = [ + [x1, y1, r1, m1], + [x2, y2, r2, m2], + [x3, y3, r3, m3], + ] + + return { + "ghosts": ghosts, + "angle_x_deg": float(angle_x * 180.0 / np.pi), + "angle_y_deg": float(angle_y * 180.0 / np.pi), + "wave_name": wave_name, + } + + def draw_on_chip(self, chip, pos_pix, r_pix, flux=1.0): + x_pix, y_pix = pos_pix.x, pos_pix.y + profile = galsim.TopHat(radius=r_pix, flux=flux) + chip.img.setOrigin(0, 0) + profile.drawImage( + image = chip.img, + center = galsim.PositionD(x_pix, y_pix), + add_to_image=True + ) + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + + +if __name__ == "__main__": + g = Ghost() + + # Original sanity checks + print(g.calculate(13, 0.0, 0.0)) + print(g.calculate(6, 0.19233, -0.24876)) + print(g.calculate(25, -0.19233, 0.24876)) diff --git a/observation_sim/mock_objects/__init__.py b/observation_sim/mock_objects/__init__.py index c6c567f9332eb3b4b109b4ff7589903b24eb7ca9..b82341546ba93c6671e3cd5748dcc53998b8356c 100755 --- a/observation_sim/mock_objects/__init__.py +++ b/observation_sim/mock_objects/__init__.py @@ -6,3 +6,4 @@ from .Star import Star from .Stamp import Stamp from .FlatLED import FlatLED from .ExtinctionMW import ExtinctionMW +from .Ghost import Ghost diff --git a/observation_sim/sim_steps/__init__.py b/observation_sim/sim_steps/__init__.py index 7669fe45bae94f3ed309840911e1d9c8660d89e7..7a48274b72aceacd430ccc26bc31b75c75106cd7 100644 --- a/observation_sim/sim_steps/__init__.py +++ b/observation_sim/sim_steps/__init__.py @@ -17,10 +17,12 @@ class SimSteps: 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, add_crosstalk from .add_LED_flat import add_LED_Flat + from .add_ghost import add_ghosts_SCI SIM_STEP_TYPES = { "scie_obs": "add_objects", + "ghost_SCI": "add_ghosts_SCI", "sky_background": "add_sky_background", "cosmic_rays": "add_cosmic_rays", "PRNU_effect": "apply_PRNU", diff --git a/observation_sim/sim_steps/add_ghost.py b/observation_sim/sim_steps/add_ghost.py new file mode 100644 index 0000000000000000000000000000000000000000..7bd4b0be61185a1b090f3e3ebcf6a0c02856c3aa --- /dev/null +++ b/observation_sim/sim_steps/add_ghost.py @@ -0,0 +1,85 @@ +import galsim + +from observation_sim.mock_objects import Ghost +from observation_sim.psf import FieldDistortion + +def add_ghosts_SCI(self, chip, filt, tel, pointing, catalog, obs_param): + # Load catalogues + if catalog is None: + self.chip_output.Log_error( + "Catalog interface class must be specified for SCIE-OBS") + raise ValueError( + "Catalog interface class must be specified for SCIE-OBS") + cat = catalog(config=self.overall_config, chip=chip, + pointing=pointing, chip_output=self.chip_output, filt=filt) + + # Apply field distortion model + if obs_param["field_dist"] is True: + fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) + else: + fd_model = None + + # Get chip WCS + if not hasattr(self, 'h_ext'): + _, _ = self.prepare_headers(chip=chip, pointing=pointing) + + chip_wcs = galsim.FitsWCS(header=self.h_ext) + ghost_model = Ghost() + + # Loop over objects + for j in range(len(cat.objs)): + obj = cat.objs[j] + + try: + sed_data = cat.load_sed(obj) + norm_filt = cat.load_norm_filt(obj) + + obj.sed, obj.param["mag_%s" % filt.filter_type.lower()], obj.param["flux_%s" % filt.filter_type.lower()] = cat.convert_sed( + mag=obj.param["mag_use_normal"], + sed=sed_data, + target_filt=filt, + norm_filt=norm_filt, + mu=obj.mu + ) + except Exception as e: + traceback.print_exc() + self.chip_output.Log_error(e) + continue + + # Select only bright stars + if obj.type != 'star' or obj.getMagFilter(filt) >= obs_param["mag_threshold"]: + continue + + # 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, ra_offset=self.ra_offset, dec_offset=self.dec_offset) + + # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane + if pos_img is None: + 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("Object missed: %s" % (obj.id)) + obj.unload_SED() + continue + + # Get number of photons for total photons for the object + nphotons_tot = obj.getElectronFluxFilt(filt, tel, pointing.exp_time) + + x_m = pos_img.x * chip.pix_size * 1e-3 + y_m = pos_img.y * chip.pix_size * 1e-3 + + ghost_list = ghost_model.calculate(chip.chipID, x_m, y_m)["ghosts"] + + for ghost in ghost_list: + x, y, r, ratio = ghost + factor = (1e-3) * chip.pix_size + x_pix = x / factor + y_pix = y / factor + r_pix = r / factor + ghost_pos = obj.getRealPos(chip.img, global_x=x_pix, global_y=y_pix, + img_real_wcs=obj.chip_wcs) + # print(f"star at {obj.getRealPos(chip.img, global_x=pos_img.x, global_y=pos_img.y, img_real_wcs=obj.chip_wcs)}") + # print("ghost info: ", ghost_pos.x, ghost_pos.y, r_pix, ratio * nphotons_tot) + ghost_model.draw_on_chip(chip, ghost_pos, r_pix, flux=ratio * nphotons_tot) + + return chip, filt, tel, pointing \ No newline at end of file