diff --git a/pyproject.toml b/pyproject.toml index 3e616be4..b505fadd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", + "romanisim", ] [project.urls] diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index 96472722..19b657b5 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -1,4 +1,4 @@ -import galsim.roman as roman +import romanisim.models as models from galsim.config import BandpassBuilder, GetAllParams, RegisterBandpassType @@ -22,10 +22,14 @@ def buildBandpass(self, config, base, logger): the constructed Bandpass object. """ req = {"name": str} - kwargs, safe = GetAllParams(config, base, req=req) + opt = {"SCA": int} + kwargs, safe = GetAllParams(config, base, req=req, opt=opt) name = kwargs["name"] - bandpass = roman.getBandpasses(red_limit=2000)[name] + if "SCA" in kwargs: + bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[name] + else: + bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] return bandpass, safe diff --git a/roman_imsim/config/hack.yaml b/roman_imsim/config/hack.yaml new file mode 100644 index 00000000..e9556791 --- /dev/null +++ b/roman_imsim/config/hack.yaml @@ -0,0 +1,175 @@ +# Default settings for roman simulation +# Includes creation of noiseless oversampled images (including PSF) +# - processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. + - roman_imsim + - galsim.roman + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. + - datetime + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: { type: ObSeqData, field: ra } + dec: { type: ObSeqData, field: dec } + pa: { type: ObSeqData, field: pa } + mjd: { type: ObSeqData, field: mjd } + + bandpass: + type: RomanBandpassTrimmed + SCA: '@image.SCA' + name: { type: ObSeqData, field: filter } + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - { type: ObSeqData, field: visit } + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 4 + mjd: { type: ObSeqData, field: mjd } + filter: { type: ObSeqData, field: filter } + exptime: { type: ObSeqData, field: exptime } + + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slightly slower than achromatic, so the difference + # is even more pronounced in this case. + draw_method: 'phot' + use_fft_bright: True + + noise: + # These are all by default turned on, but you can turn any of them off if desired: + type: RomanNoise + mjd: { type: ObSeqData, field: mjd } + stray_light: True + thermal_background: True + reciprocity_failure: True + dark_current: + turn_on: True + use_crds: True + nonlinearity: + turn_on: True + use_crds: True + ipc: + turn_on: True + use_crds: True + read_noise: + turn_on: True + use_crds: True + gain: + turn_on: True + use_crds: True + sky_subtract: False + # noise: + # type: NoNoise + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: { type: ObSeqData, field: exptime } + skip_failures: True + photon_ops: + - { type: ChargeDiff } + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + file_name: "Roman_WAS_obseq_11_1_23.fits" + visit: 12909 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: "skyCatalog/skyCatalog.yaml" + edge_pix: 512 + mjd: { type: ObSeqData, field: mjd } + exptime: { type: ObSeqData, field: exptime } + obj_types: ['diffsky_galaxy','star','snana','gaia_star'] + +output: + + nfiles: 1 + dir: output/RomanWAS_new/images/truth + file_name: + type: FormattedStr + format: "Roman_WAS_truth_%s_%i_%i.fits.gz" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + + truth: + dir: output/RomanWAS_new/truth + file_name: + type: FormattedStr + format: "Roman_WAS_index_%s_%i_%i.txt" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/roman_imsim/detector_physics.py b/roman_imsim/detector_physics.py index 7e58f8d4..ccfe4323 100644 --- a/roman_imsim/detector_physics.py +++ b/roman_imsim/detector_physics.py @@ -172,7 +172,15 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ self.pointing.sca, sky_mean=sky_mean, ) - write_fits(old_filename, new_filename, img, sky_noise, dq, self.pointing.sca, sky_mean=sky_mean) + write_fits( + old_filename, + new_filename, + img, + sky_noise, + dq, + self.pointing.sca, + sky_mean=sky_mean, + ) def get_path_name(self, use_galsim=False): @@ -309,7 +317,13 @@ def add_effects_scafile(self, im, wt, pointing): self.sky_mean, sky_noise, ) - return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise + return ( + im, + self.sky[self.sky.bounds & im.bounds] - self.sky_mean, + dq, + self.sky_mean, + sky_noise, + ) def add_effects_galsim(self, im, wt, pointing): """ @@ -439,7 +453,13 @@ def add_effects_simple(self, im, wt, pointing): self.sky_mean, sky_noise, ) - return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise + return ( + im, + self.sky[self.sky.bounds & im.bounds] - self.sky_mean, + dq, + self.sky_mean, + sky_noise, + ) def set_diff(self, im=None): if self.params["save_diff"]: @@ -509,7 +529,16 @@ def bfe(self, im): a = (a + np.fliplr(a) + np.flipud(a) + np.flip(a)) / 4.0 r = 0.5 * (3.25 / 4.25) ** (1.5) / 1.5 # source-boundary projection - B = (a[2, 2], a[3, 2], a[2, 3], a[3, 3], a[4, 2], a[2, 4], a[3, 4], a[4, 4]) + B = ( + a[2, 2], + a[3, 2], + a[2, 3], + a[3, 3], + a[4, 2], + a[2, 4], + a[3, 4], + a[4, 4], + ) A = np.array( [ @@ -584,7 +613,6 @@ def bfe(self, im): for gj in range(num_grids): for gi in range(num_grids): - a_components_pad = np.zeros( ( 4, @@ -595,7 +623,13 @@ def bfe(self, im): ) ) # (4,5,5,sub_grid,sub_grid) a_components_pad = np.zeros( - (4, 2 * nbfe + 1, 2 * nbfe + 1, bin_size * n_sub + 2 * nbfe, bin_size * m_sub + 2 * nbfe) + ( + 4, + 2 * nbfe + 1, + 2 * nbfe + 1, + bin_size * n_sub + 2 * nbfe, + bin_size * m_sub + 2 * nbfe, + ) ) # (4,5,5,sub_grid,sub_grid) for comp in range(4): @@ -858,7 +892,6 @@ def dark_current(self, im): self.im_dark = im - self.im_dark else: - dark_current_ = self.dark_current_.clip(0) # opt for numpy random geneator instead for speed @@ -989,7 +1022,6 @@ def add_persistence(self, im, pointing): im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime else: - # setup parameters for persistence Q01 = self.df["PERSIST"].read_header()["Q01"] Q02 = self.df["PERSIST"].read_header()["Q02"] @@ -1023,7 +1055,7 @@ def add_persistence(self, im, pointing): # Do linear interpolation a = np.zeros(x.shape) - a += ((x < Q01)) * x / Q01 + a += (x < Q01) * x / Q01 a += ((x >= Q01) & (x < Q02)) * (Q02 - x) / (Q02 - Q01) im.array[:, :] += a * self.df["PERSIST"][0, :, :][0] * fac_dt @@ -1049,7 +1081,7 @@ def add_persistence(self, im, pointing): a = np.zeros(x.shape) a += ((x >= Q05) & (x < Q06)) * (x - Q05) / (Q06 - Q05) - a += ((x >= Q06)) * (x / Q06) ** alpha # avoid fractional power of negative values + a += (x >= Q06) * (x / Q06) ** alpha # avoid fractional power of negative values im.array[:, :] += a * self.df["PERSIST"][5, :, :][0] * fac_dt return im @@ -1132,7 +1164,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): for dy in range(-1, 2): for dx in range(-1, 2): - array_out[ gj * grid_size : (gj + 1) * grid_size, gi * grid_size : (gi + 1) * grid_size, @@ -1258,7 +1289,6 @@ def finalize_sky_im(self, im, pointing): im = self.e_to_ADU(im) im.quantize() else: - bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096) im_pad = galsim.Image(bound_pad) im_pad.array[4:-4, 4:-4] = im.array[:, :] diff --git a/roman_imsim/noise.py b/roman_imsim/noise.py index 1d7bd711..b2fcd4ae 100644 --- a/roman_imsim/noise.py +++ b/roman_imsim/noise.py @@ -1,5 +1,5 @@ import galsim -import galsim.roman as roman +import romanisim.models as models from galsim.config import NoiseBuilder, RegisterNoiseType from astropy.time import Time @@ -38,11 +38,12 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): "stray_light": bool, "thermal_background": bool, "reciprocity_failure": bool, - "dark_current": bool, - "nonlinearity": bool, - "ipc": bool, - "read_noise": bool, + "dark_current": dict, + "nonlinearity": dict, + "ipc": dict, + "read_noise": dict, "sky_subtract": bool, + "gain": dict, } params, safe = galsim.config.GetAllParams(config, base, req={}, opt=opt, ignore=[]) @@ -51,11 +52,13 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): stray_light = params.get("stray_light", False) thermal_background = params.get("thermal_background", False) reciprocity_failure = params.get("reciprocity_failure", False) - dark_current = params.get("dark_current", False) - nonlinearity = params.get("nonlinearity", False) - ipc = params.get("ipc", False) - read_noise = params.get("read_noise", False) + # dark_current = params.get("dark_current", False) + dark_current = params.get("dark_current", {"turn_on": False, "use_crds": False}) + nonlinearity = params.get("nonlinearity", {"turn_on": False, "use_crds": False}) + ipc = params.get("ipc", {"turn_on": False, "use_crds": False}) + read_noise = params.get("read_noise", {"turn_on": False, "use_crds": False}) sky_subtract = params.get("sky_subtract", True) + gain = params.get("gain", {"turn_on": False, "use_crds": False}) base["current_noise_image"] = base["current_image"] wcs = base["wcs"] @@ -69,17 +72,17 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): # value added to sky_image. So technically, this includes things that aren't just sky. # E.g. includes dark_current and thermal backgrounds. sky_image = image.copy() - sky_level = roman.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) + sky_level = models.backgrounds.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) logger.debug("Adding sky_level = %s", sky_level) if stray_light: - logger.debug("Stray light fraction = %s", roman.stray_light_fraction) - sky_level *= 1.0 + roman.stray_light_fraction + logger.debug("Stray light fraction = %s", models.parameters.stray_light_fraction) + sky_level *= 1.0 + models.parameters.stray_light_fraction wcs.makeSkyImage(sky_image, sky_level) # The other background is the expected thermal backgrounds in this band. # These are provided in e-/pix/s, so we have to multiply by the exposure time. if thermal_background: - tb = roman.thermal_backgrounds[filter_name] * exptime + tb = models.backgrounds.thermal_backgrounds[filter_name] * exptime logger.debug("Adding thermal background: %s", tb) sky_image += tb @@ -112,39 +115,45 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): if reciprocity_failure: logger.debug("Applying reciprocity failure") - roman.addReciprocityFailure(image) + models.nonlinearity.addReciprocityFailure(img=image) - if dark_current: - dc = roman.dark_current * exptime - logger.debug("Adding dark current: %s", dc) - sky_image += dc - dark_noise = galsim.noise.DeviateNoise(galsim.random.PoissonDeviate(rng, dc)) - image.addNoise(dark_noise) + if dark_current["turn_on"]: + logger.debug("Adding dark current: %s") + # print(f'for drak_current, use_crds = {dark_current["use_crds"]}') + dc = models.DarkCurrent(usecrds=dark_current["use_crds"]) + dc.apply(img=image, exptime=exptime) - if nonlinearity: + if nonlinearity["turn_on"]: logger.debug("Applying classical nonlinearity") - roman.applyNonlinearity(image) + non_linear = models.Nonlinearity(usecrds=nonlinearity["use_crds"]) + non_linear.apply(img=image, electrons=False) # Mosby and Rauscher say there are two read noises. One happens before IPC, the other # one after. # TODO: Add read_noise1 - if ipc: + if ipc["turn_on"]: logger.debug("Applying IPC") - roman.applyIPC(image) + IPC = models.IPC(usecrds=ipc["use_crds"]) + IPC.apply(img=image) - if read_noise: - logger.debug("Adding read noise %s", roman.read_noise) - image.addNoise(galsim.GaussianNoise(rng, sigma=roman.read_noise)) + if read_noise["turn_on"]: + logger.debug("Adding read noise") + rn = models.ReadNoise(usecrds=read_noise["use_crds"]) + rn.apply(img=image) - logger.debug("Applying gain %s", roman.gain) - image /= roman.gain + if gain["turn_on"]: + logger.debug("Applying gain") + gain_model = models.Gain(usecrds=gain["use_crds"]) + gain_model.apply(img=image) # Make integer ADU now. image.quantize() if sky_subtract: logger.debug("Subtracting sky image") - sky_image /= roman.gain + if gain["turn_on"]: + logger.debug("Applying gain") + gain_model.apply(img=sky_image) sky_image.quantize() image -= sky_image diff --git a/roman_imsim/psf.py b/roman_imsim/psf.py index 033dbd10..8f6449ad 100644 --- a/roman_imsim/psf.py +++ b/roman_imsim/psf.py @@ -57,7 +57,6 @@ """ import galsim -import galsim.roman as roman from galsim.config import ( InputLoader, RegisterInputType, @@ -66,6 +65,7 @@ RegisterObjectType, ) from galsim.errors import GalSimConfigValueError +import romanisim.models as models ########################## # PSF Interpolator Input # @@ -77,6 +77,48 @@ class PSFInterpolator: """Base class for PSF interpolator""" + +class RomanPSF(object): + """Class building needed Roman PSFs.""" + + def __init__( + self, + SCA=None, + WCS=None, + n_waves=None, + bpass=None, + extra_aberrations=None, + logger=None, + ): + + logger = galsim.config.LoggerWrapper(logger) + + if n_waves == -1: + if bpass.name == "W146": + n_waves = 10 + else: + n_waves = 5 + + corners = [ + galsim.PositionD(1, 1), + galsim.PositionD(1, models.parameters.n_pix), + galsim.PositionD(models.parameters.n_pix, 1), + galsim.PositionD(models.parameters.n_pix, models.parameters.n_pix), + ] + cc = galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2) + tags = ["ll", "lu", "ul", "uu"] + self.PSF = {} + pupil_bin = 8 + self.PSF[pupil_bin] = {} + for tag, SCA_pos in tuple(zip(tags, corners)): + self.PSF[pupil_bin][tag] = self._psf_call( + SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations + ) + for pupil_bin in [4, 2, "achromatic"]: + self.PSF[pupil_bin] = self._psf_call( + SCA, bpass, cc, WCS, pupil_bin, n_waves, logger, extra_aberrations + ) + def _parse_pupil_bin(self, pupil_bin): if pupil_bin == "achromatic": return 8 @@ -86,7 +128,7 @@ def _parse_pupil_bin(self, pupil_bin): def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations): if pupil_bin == 8: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA, bpass.name, SCA_pos=SCA_pos, @@ -100,7 +142,7 @@ def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_ extra_aberrations=extra_aberrations, ) else: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA, bpass.name, SCA_pos=SCA_pos, @@ -202,10 +244,10 @@ def initPSF( self._image_xsize = image_xsize if self._image_xsize is None: - self._image_xsize = roman.n_pix + self._image_xsize = models.parameters.n_pix self._image_ysize = image_ysize if self._image_ysize is None: - self._image_ysize = roman.n_pix + self._image_ysize = models.parameters.n_pix corners = [ galsim.PositionD(1, 1), @@ -260,12 +302,12 @@ def getPSF(self, pupil_bin, pos): if pupil_bin != 8: return psf - wll = (self._image_xsize - pos.x) * (self._image_ysize - pos.y) - wlu = (self._image_xsize - pos.x) * (pos.y - 1) - wul = (pos.x - 1) * (self._image_ysize - pos.y) + wll = (models.parameters.n_pix - pos.x) * (models.parameters.n_pix - pos.y) + wlu = (models.parameters.n_pix - pos.x) * (pos.y - 1) + wul = (pos.x - 1) * (models.parameters.n_pix - pos.y) wuu = (pos.x - 1) * (pos.y - 1) return (wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"]) / ( - (self._image_xsize - 1) * (self._image_ysize - 1) + (models.parameters.n_pix - 1) * (models.parameters.n_pix - 1) ) @@ -304,6 +346,13 @@ def getKwargs(self, config, base, logger): } ignore = ["extra_aberrations"] + # If SCA is in base, then don't require it in the config file. + # (Presumably because using Roman image type, which sets it there for convenience.) + if "SCA" in base: + opt["SCA"] = int + else: + req["SCA"] = int + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) kwargs["extra_aberrations"] = galsim.config.ParseAberrations( @@ -440,7 +489,7 @@ def BuildRomanPSF(config, base, ignore, gsparams, logger): base["image_pos"], ) else: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA=base["SCA"], bandpass=base["bandpass"].name, SCA_pos=base["image_pos"], diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py index 7d362440..7524a21a 100644 --- a/roman_imsim/sca.py +++ b/roman_imsim/sca.py @@ -2,7 +2,7 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np from astropy.time import Time from galsim.config import RegisterImageType @@ -67,7 +67,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): if "bandpass" not in config: base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger) - return roman.n_pix, roman.n_pix + return models.parameters.n_pix, models.parameters.n_pix def buildImage(self, config, base, image_num, obj_num, logger): """Build an Image containing multiple objects placed at arbitrary locations. @@ -100,7 +100,7 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image.header["MJD-OBS"] = self.mjd full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat() full_image.header["FILTER"] = self.filter - full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * roman.collecting_area) + full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * models.parameters.collecting_area) base["current_image"] = full_image diff --git a/roman_imsim/scafile.py b/roman_imsim/scafile.py index de1bea8e..75574a64 100644 --- a/roman_imsim/scafile.py +++ b/roman_imsim/scafile.py @@ -4,7 +4,7 @@ import astropy.time import galsim -import galsim.roman as roman +import romanisim.models as models from galsim.config import OutputBuilder, RegisterOutputType @@ -29,7 +29,7 @@ def setup(self, config, base, file_num, logger): if "exptime" in config: base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[0] else: - base["exptime"] = roman.exptime + base["exptime"] = models.parameters.exptime # Save the detector size, so the input catalogs can use it to figure out which # objects will be visible. diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index 0b2b8c9b..19274a4e 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -3,7 +3,7 @@ """ import galsim -import galsim.roman as roman +import romanisim.models as models import numpy as np from galsim.config import ( InputLoader, @@ -68,11 +68,11 @@ def __init__( if xsize is not None: self.xsize = xsize else: - self.xsize = roman.n_pix + self.xsize = models.parameters.n_pix if ysize is not None: self.ysize = ysize else: - self.ysize = roman.n_pix + self.ysize = models.parameters.n_pix self.obj_types = obj_types self.edge_pix = edge_pix self.logger = galsim.config.LoggerWrapper(logger) @@ -184,7 +184,9 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): # Compute the flux or get the cached value. flux = ( - skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) * self.exptime * roman.collecting_area + skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) + * self.exptime + * models.parameters.collecting_area ) if np.isnan(flux): return None @@ -289,11 +291,11 @@ def SkyCatObj(config, base, ignore, gsparams, logger): message = ( "skyCatalogs selection and SCA center do not agree: \n" "skycat.sca_center: " - f"{sca_center.ra/galsim.degrees:.5f}, " - f"{sca_center.dec/galsim.degrees:.5f}\n" + f"{sca_center.ra / galsim.degrees:.5f}, " + f"{sca_center.dec / galsim.degrees:.5f}\n" "world_center: " - f"{world_center.ra/galsim.degrees:.5f}, " - f"{world_center.dec/galsim.degrees:.5f} \n" + f"{world_center.ra / galsim.degrees:.5f}, " + f"{world_center.dec / galsim.degrees:.5f} \n" f"Separation: {sep:.2e} arcsec" ) raise RuntimeError(message) diff --git a/roman_imsim/stamp.py b/roman_imsim/stamp.py index 2554d180..47fa8106 100644 --- a/roman_imsim/stamp.py +++ b/roman_imsim/stamp.py @@ -1,6 +1,6 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np from galsim.config import RegisterStampType, StampBuilder @@ -100,7 +100,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # psf = galsim.config.BuildGSObject(base, 'psf', logger=logger)[0]['achromatic'] # obj = galsim.Convolve(gal_achrom, psf).withFlux(self.flux) obj = gal_achrom.withGSParams(galsim.GSParams(stepk_minimum_hlr=20)) - image_size = obj.getGoodImageSize(roman.pixel_scale) + image_size = obj.getGoodImageSize(models.parameters.pixel_scale) # print('stamp setup3',process.memory_info().rss) base["pupil_bin"] = self.pupil_bin @@ -121,6 +121,29 @@ def setup(self, config, base, xsize, ysize, ignore, logger): return image_size, image_size, image_pos, world_pos + def buildPSF(self, config, base, gsparams, logger): + """Build the PSF object. + + For the Basic stamp type, this builds a PSF from the base['psf'] dict, if present, + else returns None. + + Parameters: + config: The configuration dict for the stamp field. + base: The base configuration dict. + gsparams: A dict of kwargs to use for a GSParams. More may be added to this + list by the galaxy object. + logger: A logger object to log progress. + + Returns: + the PSF + """ + if base.get("psf", {}).get("type", "roman_psf") != "roman_psf": + return galsim.config.BuildGSObject(base, "psf", gsparams=gsparams, logger=logger)[0] + + roman_psf = galsim.config.GetInputObj("roman_psf", config, base, "buildPSF") + psf = roman_psf.getPSF(self.pupil_bin, base["image_pos"]) + return psf + def getDrawMethod(self, config, base, logger): """Determine the draw method to use. diff --git a/roman_imsim/utils.py b/roman_imsim/utils.py index 408f8c97..19e39170 100644 --- a/roman_imsim/utils.py +++ b/roman_imsim/utils.py @@ -2,7 +2,7 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np @@ -72,9 +72,15 @@ def getPSF(self, x=None, y=None, pupil_bin=8): raise ValueError( "x,y position for pupil_bin values other than 8 not supported. Using SCA center." ) - return self.PSF.getPSF(pupil_bin, galsim.PositionD(roman.n_pix / 2, roman.n_pix / 2)) + return self.PSF.getPSF( + pupil_bin, + galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2), + ) if (x is None) | (y is None): - return self.PSF.getPSF(8, galsim.PositionD(roman.n_pix / 2, roman.n_pix / 2)) + return self.PSF.getPSF( + 8, + galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2), + ) return self.PSF.getPSF(8, galsim.PositionD(x, y)) def getWCS(self): @@ -148,7 +154,7 @@ def getPSF_Image( return psf.drawImage(self.bpass, image=stamp, wcs=wcs, method=method) photon_ops = [self.getPSF(x, y, pupil_bin)] + self.photon_ops if include_pixel and oversampling_factor > 1: - scale = galsim.roman.pixel_scale + scale = models.parameters.pixel_scale # An array (comb) of Dirac Delta functions. comb = galsim.Add( [ diff --git a/roman_imsim/wcs.py b/roman_imsim/wcs.py index 94700318..964e4d69 100644 --- a/roman_imsim/wcs.py +++ b/roman_imsim/wcs.py @@ -1,5 +1,5 @@ import galsim -import galsim.roman as roman +import romanisim.models as models from astropy.time import Time from galsim.angle import Angle from galsim.celestial import CelestialCoord @@ -7,7 +7,6 @@ class RomanWCS(WCSBuilder): - def buildWCS(self, config, base, logger): req = { @@ -21,10 +20,10 @@ def buildWCS(self, config, base, logger): kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) if "max_sun_angle" in kwargs: - roman.max_sun_angle = kwargs["max_sun_angle"] - roman.roman_wcs.max_sun_angle = kwargs["max_sun_angle"] + models.parameters.max_sun_angle = kwargs["max_sun_angle"] + models.wcs_utils.max_sun_angle = kwargs["max_sun_angle"] pointing = CelestialCoord(ra=kwargs["ra"], dec=kwargs["dec"]) - wcs = roman.getWCS( + wcs = models.wcs_utils.getWCS( world_pos=pointing, PA=kwargs["pa"], date=Time(kwargs["mjd"], format="mjd").datetime,