romanisim¶
A Roman WFI image simulator based on galsim.
We stylize the simulator Roman I-Sim and pronounce it roman - eye - sim; the package is romanisim.
Contents¶
Overview¶
romanisim simulates Roman Wide-Field Imager images of astronomical scenes described by catalogs of sources. The simulation includes:
convolution of the sources by the Roman point spread function
optical distortion
sky background
level 1 image support (3D image stack of up-the-ramp samples)
level 2 image support (2D image after calibration & ramp fitting)
point sources and analytic galaxy profiles
expected system throughput
dark current
read noise
inter-pixel capacitance
non-linearity
saturation
ramp fitting
source injection
The simulator is based on galsim and most of these features directly invoke the equivalents in the galsim.roman package. The chief additions of this package on top of the galsim.roman implementation are using “official” webbpsf PSF, and distortion and reference files from the Roman CRDS (not yet public!). This package also implements WFI up-the-ramp sampled and averaged images like those that will be downlinked from the telescope, and the official Roman WFI file format (asdf).
Future expected features include:
frame zero effects
L3 image simulations
“image-based” simulation inputs (i.e., provide an input image based on a galaxy hydro sim; romanisim applies the Roman PSF & instrumental effects on top to produce a detailed instrumental simulation)
The best way to interact with romanisim is to make an image. Running
romanisim-make-image out.asdf
will make a test image in the file out.asdf. Naturally, usually one has a
particular astronomical scene in mind, and one can’t really simulate a scene
without knowing where the telescope is pointing and when the observation is
being made. A more complete invocation would be
romanisim-make-image --catalog input.ecsv --radec 270 66 --bandpass F087 --sca 7 --date 2026 1 1 --level 1 out.asdf
where input.ecsv includes a list of sources in the scene, the
telescope boresight is pointing to (r, d) = (270, 66), the desired
bandpass is F087, the sensor is WFI07, the date is Jan 1, 2026, and a
level 1 image (3D cube of samples up the ramp) is requested.
The output of romanisim-make-image is an appropriate asdf file for
the requested level image, with the following addition. The script
adds an additional top-level branch to the asdf tree with the name
romanisim. Here’s an example:
└─romanisim (dict)
├─bandpass (str): F087
├─catalog (NoneType): None
├─date (NoneType): None
├─filename (str): out.asdf
├─level (int): 1
├─ma_table_number (int): 1
├─radec (NoneType): None
├─sca (int): 7
├─rng_seed (NoneType): None
├─simcatobj (NDArrayType): shape=(496,), dtype=void96
├─usecrds (bool): False
└─webbpsf (bool): True
These fields are simply the arguments to romanisim-make-image,
plus an additional simcatobj field with contains the x, y,
and number of photons of each simulated source.
Features not included so far:
pedestal/frame 0 features
non-linear dark features
Installation¶
To install
pip install romanisim
and you should be largely set!
There are a few dependencies that may cause more difficulty. First,
WebbPSF requires data files to
operate. See the docs
for instructions on obtaining the relevant data files and pointing the
WEBBPSF_PATH environment variable to them. This issue can be
avoided by not setting the --webbpsf argument, in which case
romanisim uses the GalSim modeling of the Roman PSF.
Second, some synthetic scene generation tools use images of galaxies
distributed separately from the main GalSim source. See here
for information on obtaining the COSMOS galaxies for use with GalSim.
The romanisim package also has a less sophisticated scene modeling
toolkit, which just renders Sersic galaxies. The command line
interface to romanisim presently uses supports Sersic galaxy
rendering, and so many users may not need to download the COSMOS galaxies.
Third, romanisim can work with the Roman CRDS system. This functionality
is not available to the general community at the time of writing.
Using CRDS requires specifying the CRDS_PATH and
CRDS_SERVER_URL variables. CRDS is not used unless the
--usecrds argument is specified; do not include this argument
unless you have access to the Roman CRDS.
That said, the basic install process looks like this:
pip install romanisim
# to get a specific version, use instead
# pip install romanisim==0.1
# to be able to run the tests for a specific version, use instead
# pip install romanisim[test]==0.1
# get webbpsf data and untar it
mkdir -p $HOME/data/webbpsf-data
cd $HOME/data/webbpsf-data
wget https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz -O webbpsf-data.tar.gz
tar -xzf webbpsf-data.tar.gz
export WEBBPSF_PATH=$PWD/webbpsf-data
# get galsim galaxy catalogs
# Note: ~5 GB each, takes a little while to download.
# Both are needed for tests. Neither are needed if you are
# exclusively using analytic model galaxies.
galsim_download_cosmos -s 23.5
galsim_download_cosmos -s 25.2
You may wish to, for example, set up a new python virtual environment before running the above, or choose a different directory for WebbPSF’s data files.
Some users report issues with the FFTW dependency of galsim on Mac Arm systems. See galsim’s installation page for hints there. In particular it may be helpful to install FFTW before galsim and romanisim.
Running the simulation¶
The primary means by which we expect most users to make images is the command line interface:
romanisim-make-image out.asdf
The combination of romanisim-make-image and various user-generated
input catalogs allows most simulator functionality to be exercised [1].
The romanisim-make-image CLI has a number of arguments to support
this functionality:
romanisim-make-image -h
usage: romanisim-make-image [-h] [--catalog CATALOG] [--radec RADEC RADEC] [--bandpass BANDPASS]
[--sca SCA] [--usecrds] [--webbpsf] [--date DATE [DATE ...]]
[--level LEVEL] [--ma_table_number MA_TABLE_NUMBER] [--seed SEED]
[--nobj NOBJ] [--boresight] [--previous PREVIOUS]
filename
Make a demo image.
positional arguments:
filename output image (fits)
optional arguments:
-h, --help show this help message and exit
--catalog CATALOG input catalog (csv) (default: None)
--radec RADEC RADEC ra and dec (deg) (default: None)
--bandpass BANDPASS bandpass to simulate (default: F087)
--sca SCA SCA to simulate (default: 7)
--usecrds Use CRDS for distortion map (default: False)
--webbpsf Use webbpsf for PSF (default: False)
--date DATE [DATE ...]
Date of observation to simulate: year month day hour minute second
microsecond (default: None)
--level LEVEL 1 or 2, for L1 or L2 output (default: 2)
--ma_table_number MA_TABLE_NUMBER
--rng_seed SEED
--nobj NOBJ
--boresight radec specifies location of boresight, not center of WFI. (default: False)
--previous PREVIOUS previous simulated file in chronological order used for persistence modeling.
(default: None)
EXAMPLE: romanisim-make-image output_image.fits
Expected arguments controlling things like the input here to simulate, the right ascension and declination of the telescope [2], the bandpass, the SCA to simulate, the level of the image to simulate (L1 or L2), the MA table to use, and the time of the observation.
Additional arguments control some details of the simulation. The
--usecrds argument indicates that reference files should be pulled
from the Roman CRDS server; this is the recommended option when CRDS
is available. The --webbpsf argument indicates that the WebbPSF package should be used to simulate
the PSF; note that this presently disables chromatic PSF rendering.
The --rng_seed argument specifies a seed to the random number
generator, enabling reproducible results.
The --nobj argument is only used when a catalog is not specified,
and controls the number of objects that are simulated in that case.
The previous argument specifies the previous simulated frame.
This information is used to support persistence
modeling.
Making images¶
Ultimately, romanisim builds image by feeding source profiles, world coordinate system objects, and point spread functions to galsim. The image and l1 modules currently implement this functionality.
The image module is responsible for translating a metadata object that specifies everything about the conditions of the observation into objects that the simulation can understand. The metadata object follows the metadata that real WFI images will include; see here for more information.
The parsed metadata is used to make a counts image that is an idealized image containing the number of photons each WFI pixel would collect over an observation. It includes no systematic effects or noise beyond Poisson noise from the astronomical scene and backgrounds. Actual WFI observations are more complicated than just noisy versions of this idealized image, however, for several reasons:
WFI pixels have a very uncertain pedestal.
WFI pixels are sampled “up the ramp” during an observation, so a number of reads contribute to the final estimate for the rate of photons entering each pixel.
WFI reads are averaged on the telescope into resultants; ground images see only resultants.
These idealized count images are then used to either make a level 2 image or a level 1 image, which are intended to include the effects of these complications. The construction of L1 images is described here, and the construction of L2 images is described here.
romanisim.image Module¶
Roman WFI simulator tool.
Based on galsim’s implementation of Roman image simulation. Uses galsim Roman modules for most of the real work.
Functions¶
|
Add sources to an image. |
|
Gather reference data corresponding to metadata. |
|
Filter sources to those landing on an image. |
|
Wrap a galsim simulated image with ASDF/roman_datamodel metadata. |
|
Simulate an image in a filter given resultants. |
|
This routine kicks the tires on everything in this module. |
|
Simulate a sequence of observations on a field in different bandpasses. |
|
Simulate total counts in a single SCA. |
|
Add some simulated counts to an image. |
|
Trim a Table of objects down to those falling near an image. |
Making L1 images¶
An L1 (level 1) image is a “raw” image received from the detectors. The actual measurements made on the spacecraft consist of a number of non-destructive reads of the pixels of the H4RG detectors. These reads have independent read noise but because the pixels count the total number of photons having entered each pixel, the Poisson noise in different reads of the same pixel is correlated.
Because the telescope has limited bandwidth, every read is not transferred to ground stations. Instead, reads are averaged into “resultants” according to a specification called a MultiAccum table, and these resultants are transferred, archived, and analyzed. These resultants make up an L1 image, which romanisim simulates.
L1 images are created using an idealized counts image described here, which contains the number of photons each pixel of the detector would receive absent any instrumental systematics. To transform this into an L1 image, these counts must be apportioned into reads and averaged into resultants, and instrumental effects must be added.
This process proceeds by simulating each read, drawing the appropriate number of photons from the total number of photons for each read following a binomial distribution. These photons are added to a running sum that is then averaged into a resultant according to the MultiAccum table specification. This process requires drawing random numbers from the binomial distribution for every read of every pixel, and so can take on the order of a minute, but it allows detailed simulation of the statistics of the noise in each resultant together with their correlations. It also makes it straightforward to add various instrumental effects into the simulation accurately, since these usually apply to individual reads rather than to resultants (e.g., cosmic rays affect individual reads, and their affect on a resultant depends on the read in the resultant to which they apply).
After apportioning counts to resultants, systematic effects are added to the resultants. Presently only read noise is added. The read noise is averaged down like \(1/\sqrt{N}\), where \(N\) is the number of reads contributing to the resultant.
Nonlinearity¶
Non-linearity is considered when L1 images are constructed and a non-linearity model is provided (e.g., from CRDS). We treat non-linearity as a difference between the electrons captured in the detector and the amount of signal read out. This function is modeled as a high order polynomial, and the coefficients of this polynomial and its inverse are stored in CRDS (linearity, inverselinearity reference files for each detector). When assigning counts to each read, these are transformed through the inverselinearity polynomial for each pixel and then added to the resultant buffer to account for this effect. The linearity polynomial then corrects for this effects as part of calibrating an L1 file and eventually producing an L2 file.
The linearity polynomials are occasionally poorly determined or cannot be computed. When marked as problematic in the reference file, we use trivial polynomials (i.e., the identity), and mark the affected pixels with a DQ bit indicating a problematic linearity correction.
Interpixel Capacitance¶
Interpixel capacitance (IPC) is added following non-linearity and before read-out. Read noise remains independent among different pixels but the Poisson noise is correlated between pixels by the IPC. We simply convolve the resultants by a 3x3 kernel after apportioning counts to resultants and applying non-linearity but before adding read noise.
This is slightly different than including IPC in the PSF kernel because including IPC in the PSF kernel leaves the Poisson noise uncorrelated.
Persistence¶
Persistence is implemented in the simulator following Sanchez+2023. This follows the Fermi description of persistence implemented in GalSim, where the flux in electrons per second recorded in a pixel is parameterized in terms of the total number of counts recorded in an earlier frame.
Here \(P(x, t)\) is the rate in electrons per second that the pixel records \(t\) seconds following receiving a total number of electrons \(x\). The parameters \(A\), \(x_0\), \(\delta x\), \(\alpha\), \(\gamma\) may vary from pixel to pixel, though are presently fixed to global constants. This equation for the rate only applies to pixels which were illuminated more than to fill more than their half-well. We follow GalSim and linearly increase the persistence from 0 to the half-well value for illuminations between 0 and half-well.
This persistence rate is sampled with a Poisson distribution and added to each pixel read-by-read and incorporated into the resultants in the L1 images.
Persistence-affected pixels are expected to be rare, and are tracked sparsely via a list of the indices of affected pixels, the amount of the illumination, and the times of their illumination. Pixels are dropped from persistence tracking when their persistence rate is less than one electron per 100 seconds. If the same pixel is receives large fluxes multiple times, these are treated as two independent events and the resulting persistence flux is handled by summing the persistence rates given above over each event.
Cosmic rays¶
Cosmic rays are added to the simulation read-by-read. The cosmic ray parameters follow Wu et al. (2023). The locations of cosmic rays are chosen at random to sample the focal plane uniformly. Lengths are chosen according to a power law distribution \(p(l) \\sim l^{-4.33}\), with lengths between 10 and 10,000 microns. Charge deposition rates per micron are selected from a Moyal distribution located at 120 electrons per micron with a width of 50 electrons per micron. An idealized charge is computed for each pixel in a read according to the product of the deposition rate per micron and the length of the cosmic ray’s path within that pixel. This idealized charge is Poisson sampled and added to the relevant pixels in a read.
romanisim.l1 Module¶
Convert images into L1 images, with ramps.
We imagine starting with an image that gives the total number of counts from all Poisson processes (at least: sky, sources, dark current). We then need to redistribute these counts over the resultants of an L1 image.
The easiest thing to do, and probably where I should start, is to sample the image read-by-read with a binomial draw from the total counts weighted by the chance that each count landed in this particular time window. Then gather those for each resultant and average, as done on the spacecraft.
It’s tempting to go straight to making the appropriate resultants. Following Casertano (2022?), the variance in each resultant is:
where f is the count rate, N is the number of reads in the resultant, and \(\tau\) is the ‘variance-based resultant time’
where the t_k is the time of the kth read in the resultant.
For uniformly spaced reads,
where t_0 is the time of the first read in the resultant and d is the spacing of the reads.
So that gives the variance from Poisson counts in resultant. But how should we draw random numbers to get that variance and the right mean? I can separately control the mean and variance by scaling the Poisson distribution, but I’m not sure that’s doing the right thing with the higher order moments.
It probably isn’t that expensive to just work out all of the reads, individually, which will also allow more natural incorporation of cosmic rays down the road. So let’s take that approach instead for the moment.
How do we want to specify an L1 image? An L1 image is defined by a total count image and a list of lists \(t_{i, j}\), where \(t_{i, j}\) is the time at which the jth read in the ith resultant is made. We demand \(t_{i, j} > t_{k, l}\) whenever i > k or whenever i = k and j > l.
Things this doesn’t allow neatly:
jitter in telescope pointing: the rate image is the same for each read/resultant
weird non-linear systematics in darks?
Some systematics need to be applied to the individual reads, rather than to the final image. Currently linearity, persistenc, and CRs are implemented at individual read time. I need to think about when in the chain things like IPC, etc., come in. But it still seems correct to first generate the total number of counts that an ideal detector would measure from sources, and then apply these effects read-by-read as we build up the resultants—i.e., I expect the current framework will be able to handle this without major issue.
This approach is not super fast. For a high latitude set of resultants, generating all of the random numbers to determine the apportionment takes 43 s on the machine I’m currently using; this will scale linearly with the number of reads. That’s longer than the actual image production for the dummy scene I’m using (though only ~2x longer).
I don’t have a good way to speed this up. Explicitly doing the Poisson noise on each read from a rate image (rather than the apportionment of an image that already has Poisson noise) is 2x slower—generating billions of random numbers just takes a little while.
Plausibly I could figure out how to draw numbers directly from what a resultant is rather than looking at each read individually. That would likely bring a ~10x speed-up. The read noise there is easy. The poisson noise is a sum of scaled Poisson variables:
where \(c_i\) is a Poisson-distributed variable. The sum of Poisson-distributed variables is Poisson-distributed, but I wasn’t immediately able to find anything about the sum of scaled Poisson-distributed variables. The result is clearly not Poisson-distributed, but maybe there’s some special way to sample from that directly.
If we did sample from that directly, we’d still also need to get the sum of the counts in the reads comprising the resultant. So I think you’d need a separate draw for that, conditional on the number you got for the resultant. Or, reversing that, you might want to draw the total number of counts first, e.g., via the binomial distribution, and then you’d want to draw a number for what the average number of counts was among the reads comprising the resultant, conditional on the total number of counts. Then
is some kind of statistic of the multinomial distribution. That sounds a little more tractable?
We want to draw from \(\sum (N-i) c_i\). I think the probabilities are always \(1/N\), with the possible small but important exception of ‘skipped’ or ‘dropped’ reads, in which case the first read would be more like \(2/(N+1)\) and all the others \(1/(N+1)\). If the probabilities were always 1/N, this looks vaguely like it could have a nice analytic solution. Otherwise, I don’t immediately see a route forward. So I am not going to pursue this avenue further.
Functions¶
|
Add IPC to resultants. |
|
Adds read noise to resultants. |
|
Apportion counts to resultants given read times. |
|
Package and optionally write out an L1 frame. |
|
Make an L1 image from a counts image. |
|
Get the times of each read going into resultants for a read_pattern. |
|
Convert a set of times tij to corresponding probabilities for sampling. |
|
Verify that a set of times tij for a valid resultant. |
romanisim.nonlinearity Module¶
Routines to handle non-linearity in simulating ramps.
The approach taken here is straightforward. The detector is accumulating photons, but the capacitance of the pixel varies with flux level and so the mapping between accumulated photons and read-out digital numbers changes with flux level. The CRDS linearity and inverse-linearity reference files describe the mapping between linear DN and observed DN. This module implements that mapping. When simulating an image, the photons entering each pixel are simulated, and then before being “read out” into a buffer, are transformed with this mapping into observed counts. These are then averaged and emitted as resultants.
Functions¶
|
Correct the observed counts for non-linearity. |
|
Fix cases of zeros and NaNs in non-linearity coefficients. |
Classes¶
|
Keep track of non-linearity and inverse non-linearity coefficients. |
Class Inheritance Diagram¶
romanisim.persistence Module¶
Persistence module.
This module implements a persistence simulation following Sanchez+2023.
Functions¶
|
The Fermi model for persistence: A * (x/x0)**alpha * (t/1000.)**(-gamma) / (exp(-(x-x0)/dx) + 1) For influence level below the half well, the persistence is linear in x. |
Classes¶
|
Track persistence information. |
Class Inheritance Diagram¶
romanisim.cr Module¶
Functions¶
|
A function for performing inverse transform sampling. |
|
Return unnormalized Moyal distribution, which approximates a Landau distribution and is used to describe the energy loss probability distribution of a charged particle through a detector. |
|
Return unnormalized power-law distribution parameterized by a log-log slope, used to describe the cosmic ray path lengths. |
|
Generates cosmic ray parameters randomly sampled from distribution. |
|
Adds CRs to an existing image. |
|
Given a starting and ending pixel, returns a list of pixel coordinates (ii, jj) and their traversed path lengths. |
L2 images¶
L2 images are constructed from L1 images, which are in turn constructed from idealized count images. This means that even when constructing L2 images, one must go through the process of simulating how counts get apportioned among the various reads. It is challenging to realistically model the statistics in the noise of L2 images without going through this process.
L2 images are constructed by doing “ramp fitting” on the level 1 images. Each pixel of a level 1 image is a series of “resultants”, giving the measured value of that pixel averaged over a series on non-destructive reads as the exposure is being observed. A simple model for a pixel is that its flux as a function of time is simply two numbers: a pedestal and a linear ramp representing the rate at which photons are detected by the pixel. Ramp fitting turns these one-dimensional series of resultants into a “slope” image that is of interest astronomically. Due to details of the H4RG detectors, the pedestals of the ramp vary widely from exposure to exposure, and so current fitting completely throws away as non-astronomical any information in the ramp that is sensitive to the pedestal.
The package contains two algorithms for ramp fitting. The first uses “optimal” weighting, considering the full covariance matrix between each of the resultants stemming from read & Poisson noise from the sources. The covariance is inverted and combined with the design matrix in the usual least-squares approach to solve for the optimal slope and pedestal measurements for each pixel. This approach is naively expensive, but because the covariance matrices for each pixel for a one-dimensional family depending only on the ratio of the flux in the pixel to the read variance, the relevant matrices can be precomputed. These are then interpolated between for each pixel and summed over to get the parameters and variances for each pixel.
This approach does not handle cosmic rays or saturated pixels well, though for modest sized sets of resultants introducing an additional series of fits for the roughly \(2 n_\mathrm{resultant}\) sub-ramps would be straightforward. That approach would also naturally handle saturated ramps. Even explicitly computing every possible \(n_\mathrm{resultant} (n_\mathrm{resultant} - 1) / 2\) subramp would likely still be quite inexpensive for modestly sized ramps.
The second approach follows Casertano+2022. In this approach, a diagonal set of weights is used in place of the full covariance matrix. The choice of weights depend on the particular pattern of reads assigned to each resultant and the amount of flux in the ramp, allowing them to interpolate from simply differencing the first and last resultants when the flux is very large to weighting the resultants by the number of reads when the flux is zero. This approach more efficiently handles dealing with ramps that have been split by cosmic rays, and obtaining uncertainties within a few percent of the “optimal” weighting approach. For these cases, we report final ramp slopes and variances derived from the inverse variance weighted subramp slopes and variances, using the read-noise derived variances.
This is a fairly faithful representation of how level two image construction works, so there are not many additional effects to add here. But mentioning some limitations:
We have a simplistic saturation treatment, just clipping resultants that exceed the saturation level to the saturation level and throwing a flag.
romanisim.ramp Module¶
Ramp fitting routines.
The simulator need not actually fit any ramps, but we would like to do a good job simulating the noise induced by ramp fitting. That requires computing the covariance matrix coming out of ramp fitting. But that’s actually a big part of the work of ramp fitting.
There are a few different proposed ramp fitting algorithms, differing in their weights. The final derived covariances are all somewhat similarly difficult to compute, however, since we ultimately end up needing to compute
for the “optimal” case, or
for some alternative weighting.
We start trying the “optimal” case below.
For the “optimal” case, a challenge is that we don’t want to compute \(C^{-1}\) for every pixel individually. Fortunately, we only need \((A^T C^{-1} A)^{-1}\) (which is only a 2x2 matrix) for variances, and only \((A^T C^{-1} A)^{-1} A^T C^{-1}\) for ramp fitting, which is 2xn. Both of these matrices are effectively single parameter families, depending after rescaling by the read noise only on the ratio of the read noise and flux.
So the routines in these packages construct these different matrices, store them, and interpolate between them for different different fluxes and ratios.
Functions¶
|
Constructs covariance matrix for first finite differences of unevenly sampled resultants. |
|
Construct the \(k_i\) weights and variances for ramp fitting. |
|
Construct \(A^T C^{-1} A\) and \(A^T C^{-1}\), the matrices needed to fit ramps from resultants. |
|
Fit ramps following Casertano+2022, including averaging partial ramps. |
|
Fit ramps following Casertano+2022, only using full ramps. |
|
Construct a grid of \(k\) and covariances for the values of flux_on_readvar. |
|
Construct the tau for each resultant from a read_pattern. |
|
Construct the mean times for each resultant from a read_pattern. |
|
Convert resultants to their finite differences. |
|
Simulate many ramps with a particular flux, read noise, and read_pattern. |
Classes¶
|
Ramp fitting tool aiding efficient fitting of large number of ramps. |
Class Inheritance Diagram¶
Reference files¶
romanisim uses reference files from CRDS in order to simulate realistic images. The following kinds of reference files are used:
read noise
dark current
flat field
gain
distortion map
The usage of these is mostly straightforward, but we provide here a few notes.
Read Noise¶
The random noise on reading a sample contributing to a ramp in an L1 image is scaled by the read noise reference file.
Dark Current¶
CRDS provides dark current images for each possible MA table, including the averaging of the dark current into resultants. This simplifies subtraction from L1 images and allows effects beyond a simple Poisson sampling of dark current electrons in each read. But it’s unwieldy for a simulator because any effects beyond simple Poisson sampling of dark current electrons are not presently defined well enough to allow simulation. So the simulator simply takes the last resultant in the dark current resultant image and scales it by the effective exposure time of that resultant to get a dark current rate. This rate then goes into the idealized “counts” image which is then apportioned into the reads making up the resultants of an L1 image.
Flat field¶
Implementation of the flat field requires a little care due to the desire to support galsim’s “photon shooting” rendering mode. This mode does not create noise-free images but instead only simulates the number of photons that would be actually detected in a device. We want to start by simulating the number of photons each pixel would record for a flat field of 1, and then sample that down by a fraction corresponding to the actual QE of each pixel. That works fine supposing that the flat field is less than 1, but does not work for values of the flat field greater than 1. So we instead do the initial galsim simulations for a larger exposure time than necessary, scaled by the maximum value of the flat field, and then sample down by \(\mathrm{flat}/\mathrm{maxflat}\). That’s all well and good as long as there aren’t any spurious very large values in the flat field. I haven’t actually seen any such values yet and so haven’t tried to address that case (e.g., by clipping them).
Gain¶
Photons from the idealized “counts” image are scaled down to ADU before quantization during L1 creation, and then converted back to electrons before ramp fitting when making L2 images.
Distortion map¶
World coordinate systems for WFI images are created by placing the telescope boresight at \(\mathrm{V2} = \mathrm{V3} = 0\), and then applying the distortion maps from CRDS to convert from V2V3 to pixels.
Catalogs¶
The simulator takes catalogs describing objects in a scene and generates images of that scene. These catalogs have the following form:
ra dec type n half_light_radius pa ba F087
float64 float64 str3 float64 float64 float64 float64 float64
------- ------- ---- ------- ----------------- ------- ------- --------
269.9 66.0 SER 1.6 0.6 165.6 0.9 1.80e-09
270.1 66.0 SER 3.6 0.4 71.5 0.7 3.35e-09
269.8 66.0 PSF -1.0 0.0 0.0 1.0 2.97e-10
269.9 66.0 SER 2.5 0.8 308.8 0.7 1.50e-09
269.8 65.9 SER 3.9 0.9 210.0 0.9 3.28e-10
270.1 66.0 SER 4.0 1.1 225.1 1.0 1.61e-09
269.9 65.9 SER 1.5 0.3 271.8 0.6 1.13e-09
269.9 65.9 SER 2.9 2.3 27.6 1.0 3.28e-09
269.9 66.0 SER 1.1 0.3 4.3 1.0 9.99e-10
The following fields must be specified for each source:
ra: the right ascension of the source
dec: the declination of the source
type: PSF or SER; whether the source is a point source or Sersic galaxy
n: the Sersic index of the source. This value is ignored if the source is a point source.
half_light_radius: the half light radius of the source in arcseconds. This value is ignored if the source is a point source.
pa: the position angle of the source, in degrees east of north. This value is ignored if the source is a point source.
ba: the major-to-minor axis ratio. This value is ignored if the source is a point source.
Following these required fields is a series of columns giving the fluxes of the the sources in “maggies”; the AB magnitude of the source is given by \(-2.5*\log_{10}(\mathrm{flux})\). In order to simulate a scene in a given bandpass, a column with the name of that bandpass must be present giving the total fluxes of the sources. Many flux columns may be present, and other columns may also be present but will be ignored.
The simulator then renders these images in the scene and produces the simulated L1 or L2 images.
The simulator API includes a few simple tools to generate parametric distributions of stars and galaxies. The make_stars and make_galaxies routines make random catalogs of stars and galaxies. The number of stars and galaxies can be adjusted. Likewise, the power law index by which the sources’ magnitudes are sampled can be adjusted, as can their limiting magnitudes. Galaxy Sersic parameters, half-light radii, and position angles are chosen at random, with a rough attempt to make brighter galaxies appropriately larger (i.e., conserving surface brightness). Stars can be chosen to be distributed with a King profile. This functionality is however very rudimentary and limited, and is better suited for toy problems than real scientific work. We expect scientific uses to be driven by custom-created catalogs rather than these simple routines.
romanisim.catalog Module¶
Catalog generation and reading routines.
This module provides basic routines to allow romanisim to render scenes based on catalogs of sources in those scenes.
Functions¶
|
Make a dummy catalog for testing purposes. |
|
Make a dummy table catalog. |
|
Make a simple parametric catalog of galaxies. |
|
Make a simple parametric catalog of stars. |
|
Read a catalog into a list of CatalogObjects. |
|
Read a astropy Table into a list of CatalogObjects. |
Classes¶
|
Simple class to hold galsim positions and profiles of objects. |
Class Inheritance Diagram¶
APT file support¶
The simulator possesses rudimentary support for simulating images from APT files. In order to simulate a scene, romanisim needs to know what’s in the scene, as specified by a catalog. It also needs to know where the telescope is pointed, the roll angle of the telescope, the date of the observation, and the bandpass. Finally, it needs to know what the MultiAccum table of the observation is—roughly, how long the exposure is and how the reads of the detector should be averaged into resultants.
Much of this information is available in an APT file. A rudimentary APT file reader can pull out the right ascension and declinations of observations, as well as the filters requested. However, support for roll angles is not yet included. APT files do not include the dates of observation, so this likewise is not included. APT files naturally do not contain catalogs of sources in the field, so some provision must be made for adding this information.
This module is not yet fully baked.
romanisim.apt Module¶
Very simple APT reader.
Converts an APT file into a list of (ra, dec, angle, filter, date, exposure time) needed for generating observations. This is adequate for reading in a few of the example Roman APTs but only supports a tiny fraction of what an APT file seems able to do.
Functions¶
|
Read an APT file, returning a list of observations. |
Classes¶
|
An observation of a target. |
|
A target for observation. |
Class Inheritance Diagram¶
Bandpasses¶
The simulator can render scenes in a number of different bandpasses. The choice of bandpass affects the point spread function used, the sky backgrounds, the fluxes of sources, and the reference files requested.
At present, romanisim simply passes the choice of bandpass to other packages—to webbpsf for PSF modeling, to galsim.roman for sky background estimation, to CRDS for reference file selection, or to the catalog for the selection of appropriate fluxes. However, because catalog fluxes are specified in “maggies” (i.e., in linear fluxes on the AB scale), the simulator needs to know how to convert between a maggie and the number of photons Roman receives from a source. Accordingly, the simulator knows about the AB zero points of the Roman filters, as derived from https://roman.gsfc.nasa.gov/science/WFI_technical.html .
One technical note: it is unclear what aperture is used for the bandpasses provided by Goddard. The Roman PSF formally extends to infinity and some light is received by the detector but is so far from the center of the PSF that it is not useful for flux, position, or shape measurements. Often for the purposes of computing effective area curves only light landing within a fixed aperture is counted. Presently the simulator assumes that an infinite aperture is used. This can result in an approximately 10% different flux scale than more reasonable aperture selections.
romanisim.bandpass Module¶
Roman bandpass routines
The primary purpose of this module is to provide the number of counts per second expected for sources observed by Roman given a source with the nominal flat AB spectrum of 3631 Jy. The ultimate source of this information is https://roman.gsfc.nasa.gov/science/WFI_technical.html .
Functions¶
|
Compute the AB zero point fluxes for each filter. |
|
Compute the count rate in a given filter, for a specified SED. |
|
Get the zero point flux for a particular bandpass. |
|
Read an effective area file from Roman. |
Point Spread Function Modeling¶
The simulator has two mechanisms for point source modeling. The first uses the galsim implementation of the Roman point spread function; for more information, see the galsim Roman documentation. The second uses the webbpsf package to make a model of the Roman PSF.
In the current implementation, the simulator uses a linearly varying, achromatic bandpass for each filter when using webbpsf. That is, the PSF does not vary depending on the spectrum of the source being rendered. However, it seems straightforward to implement either of these modes in the context of galsim, albeit at some computational expense.
When using the galsim PSF, galsim’s “photon shooting” mode is used for efficient rendering of chromatic sources. When using webbpsf, FFTs are used to to do the convolution of the intrinsic source profile with the PSF and pixel grid of the instrument.
romanisim.psf Module¶
Roman PSF interface for galsim.
galsim.roman has an implementation of Roman’s PSF based on the aperture and some estimates for the wavefront errors over the aperture as described by amplitudes of various Zernicke modes. This seems like a very good approach, but we want to add here a mode using the official PSFs coming out of webbpsf, which takes a very similar overall approach.
galsim’s InterpolatedImage class makes this straightforward. Future work should consider the following:
how do we want to deal with the dependence of the PSF on the source SED? It’s possible I can just subclass ChromaticObject and implement evaluateAtWavelength, possibly also stealing the _shoot code from ChromaticOpticalPSF?
Functions¶
|
Make a PSF profile for Roman at a specific detector location. |
|
Make a PSF profile for Roman. |
Classes¶
|
Spatially variable PSF wrapping GalSim profiles. |
Class Inheritance Diagram¶
World Coordinate Systems & Distortion¶
The simulator has two options for modeling distortion and world coordinate systems. The first is to use the routines in the galsim.roman package; see GalSim’s documentation for more information. The second is to use distortion reference files from the Calibration References Data System (CRDS).
The latter system works by grabbing reference distortion maps for the appropriate detector and filter combinations from the Roman CRDS server. These distortion maps are then wrapped into a galsim WCS object and fed to galsim’s rendering pipeline.
romanisim.wcs Module¶
Roman WCS interface for galsim.
galsim.roman has an implementation of Roman’s WCS based on some SIP coefficients for each SCA. This is presumably plenty good, but here we take the alternative approach of using the distortion functions provided in CRDS. These naturally are handled by the gWCS library, but galsim only naturally supports astropy WCSes via the ~legacy interface. So this module primarily makes the wrapper that interfaces gWCS and galsim.CelestialWCS together.
This presently gives rather different world coordinates given a specific telescope boresight. Partially this is not doing the same roll_ref determination that galsim.roman does, which could be fixed. But additionally the center of the SCA looks to be in a different place relative to the boresight for galsim.roman than for what I get from CRDS. This bears more investigation.
Functions¶
|
Convert a GalSim WCS object into a GWCS object. |
|
Add WCS info to parameters dictionary. |
|
Get a WCS object for a given sca or set of CRDS parameters. |
|
Create a gWCS from a target position, a roll, and a distortion map. |
|
Convert a FITS WCS to a GWCS. |
Classes¶
|
This WCS uses gWCS to implent a galsim CelestialWCS. |
Class Inheritance Diagram¶
Parameters¶
The parameters module contains useful constants describing the Roman telescope. These include:
the default read noise when CRDS is not used
the number of border pixels
the specification of the read indices going into resultants for a fiducial L1 image for a handful of MA tables
the read time
the default saturation limit
the default IPC kernel
the definition of the reference V2/V3 location in the focal plane to which to place the given ra/dec
the default persistence parameters
the default cosmic ray parameters
These values can be overridden by specifying a yaml config file on the command line to romanisim-make-image.
romanisim.parameters Module¶
Parameters class storing a few useful constants for Roman simulations.
Utilities¶
The simulator utility module consists of miscellaneous utility routines intended to be of broad use. Present examples include:
turning astropy coordinates into galsim coordinates and vice-versa
making an RGB image from image slices
generating points at random in a region on the sky
romanisim.util Module¶
Miscellaneous utility routines.
Functions¶
|
Fill out the metadata dictionary, modifying it in place. |
|
Turn a SkyCoord into a CelestialCoord. |
|
Compute the King (1962) profile. |
|
Choose locations at random at given radii from coord. |
|
Choose locations at random within radius of coord. |
|
Sample points from a King distribution |
|
Sample distances from a King (1962) profile. |
|
Scales three flux images into a range of luminosity for displaying. |
|
Turn a CelestialCoord into a SkyCoord. |